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% : Abstract 

We present a lattice method to compute bubble nucleation rates at radiatively 

induced first order phase transitions, in high temperature, weakly coupled 

>^ • field theories, nonperturbatively. A generalization of Langer's approach, it 

CN . makes no recourse to saddle point expansions and includes completely the 

I , dynamical prefactor. We test the technique by applying it to the electroweak 

0^ \ phase transition in the minimal standard model, at an unphysically small 

3^ ' Higgs mass which gives a reasonably strong phase transition {X/g"^ = 0.036, 

O ■ which corresponds to mn /rny/ = 0.54 at tree level but does not correspond 

r~| ' to a positive physical Higgs mass when radiative effects of the top quark are 

Mh, included), and compare the results to older perturbative and other estimates. 

Qh! While two loop perturbation theory slightly under-estimates the strength of 

the transition measured by the latent heat, it over-estimates the amount of 

supercooling by a factor of 2. 






X 

^ , I. INTRODUCTION 

Electroweak baryogenesis provides one of the best motivated and most testable mecha- 
nisms for the origin of the cosmological baryon number abundance. However, several ingredi- 
ents are missing before we can make quantitative predictions. One set of needed ingredients 
is particle physics inputs. For instance, it is difficult to say much about electroweak matter 
at temperatures T ~ lOOGeV, where electroweak symmetry is "restored" [| and baryogenesis 
may occur, when we still do not know the Higgs mass. It is also absolutely necessary to 
know what other light (m <, 150GeV) scalars there are, what couplings they have with the 



^There is no qualitative distinction between the "symmetric" and "broken" electroweak phases, 
which are in fact analytically connected, and the symmetry is never truly "broken" even in vacuum; 
but when the phase transition is reasonably strong their quantitative behavior is very different. 
We will use the "symmetric" and "broken" terminology because it is convenient and widespread. 



Higgs boson(s), and what CP violation is operative at electroweak energy and temperature 
scales. Answering these questions will require new experimental results, and we will not 
have more to say about that here. 

However, even if we knew the complete electroweak theory we could not at this time make 
accurate predictions of what baryon number would be cosmologically produced, because we 
do not have a complete set of reliable computational tools for studying the electroweak phase 
transition and the physical processes responsible for baryogenesis. This is best illustrated 
by briefly reviewing what the scenario is, and which aspects we do or do not currently have 
good control of. 

Assuming a standard thermal history for the early universe back to0 T ~ lOOGeV, 
electroweak baryogenesis appears to be possible only if there is a fairly strong first order 
electroweak phase transition 0,^. The electroweak phase transition, if there is one, is 
radiatively induced, and determining its order and strength is a difficult problem with a 
long history. In summary, perturbation theory proves to be of some limited use when the 
phase transition is strong, but a reliable calculation of the strength of the electroweak phase 
transition requires a nonperturbative lattice calculation. The equipment for performing an 
accurate lattice calculation now exists, using either a 3-dimensional effective theory ||^-§| or 
4-dimensional SU(2) gauge + Higgs theory P,pD|, so this part of the problem is solved. 



If the electroweak phase transition is first order, then the universe will remain in the 
"symmetric" phase even after it is no longer thermodynamically favored. How deeply it 
supercools is the topic of this paper and we will return to it momentarily. After sufficiently 
deep supercooling, critical bubbles of the broken phase form at a cosmologically relevant rate, 
expand, and coalesce, completing the phase transition.^ It is the expansion of these bubbles 
into the symmetric phase which is expected to generate the baryon number. Specifically, 
the moving phase interface (bubble wall) can inject a CP violating fiux of particles into 
the symmetric phase [|12|, where baryon number violation is efficient |ll3| , |l4| . Recently, the 



efficiency of the baryon number violation has been pinned down fairly accurately |T^-|Ty[. 
However, neither the expansion of bubbles into the symmetric phase, nor the generation and 
propagation of CP violating particle fluxes, can yet be calculated with much precision or 
confldence, though there has been recent progress on both problems [pO|-|23[ . 



The bubble nucleation rate enters the flnal baryon number asymmetry by determining 
the amount of supercooling which occurs. The more supercooling, the greater the free energy 
difference between the phases, and the faster phase interfaces propagate. This in turn would 
mean a larger injected CP violating flux (since the flux must vanish in equilibrium in a 
CPT respecting theory). However, it would also mean that the phase interface would more 
quickly catch back up with particles it injected into the symmetric phase, which could have 
the reverse effect. The detailed dependence on supercooling and bubble wall velocity may 
be complicated, see for instance [p4| , ^ . How deep the supercooling proceeds can also be 



^See Q for an interesting discussion of what could happen if this assumption were not true. 

^It has been argued that the phase transition can also proceed by coalescence of "subcritical 
bubbles" with almost no supercooling [11|. We feel our technique and results, presented here, put 



this idea to rest for the phase transition strength we are interested in. 



important because the universe heats during the phase transition, as the latent heat of the 
symmetric phase is released. We find below that the supercooling is less than in perturbation 
theory, while the latent heat is sometimes more; so this effect is more important than one 
might have anticipated. For the parameters we will study, there is enough latent heat, and 
little enough supercooling, that the universe reheats to the equilibrium temperature Teq 
before all space has converted to the broken phase. (The remaining space would convert 
much more slowly, as the expansion of the universe continues to absorb heat from the 
plasma.) 

We are interested in a regime where the bubble nucleation rate is extremely small. This 
is because the phase transition completes when the bubble nucleation rate is, very roughly, 
around one bubble per Hubble volume per Hubble time. But at T ~ lOOGeV, a Hubble 
time is thubbie ~ H~^ ~ mpi/T^ ~ 10^^/T; so the rate of bubble nucleations must be 
~ (10~^^T)^ = IQ-^^T^ = e~^^^T^. A more careful calculation, accounting for how much 
time the phase transition takes to complete, shows that the nucleation rate must be about 
g-i062i4 '\Yhen the nucleation rate is so small, it is "almost" a thermodynamical quantity, 
set by the free energy of the "critical bubble." This gives us a hint at how to determine it 
on the lattice; we must determine the free energy of a critical bubble. However, deciding 
what precisely this means and how to go about doing it, and relating the result to the real 
time rate for a rare process, require some work. In this paper we present a quantitative 
approach to address this problem, and we carry out our program for the minimal standard 
model, at zero Weinberg angle and an unphysical Higgs mass. Clearly this case will not be 
of direct physical interest. However, it allows us to determine how well the technique works, 
and to compare the bubble nucleation rate to what we would get using one of several less 
rigorous methods, such as a perturbative calculation of the critical bubble free energy or a 
"thin wall" treatment from either perturbative or nonperturbative inputs. We find that the 
"thin wall" treatment gives a reasonable answer but is not extremely accurate, while the 
perturbative approach is quite bad unless Higgs field wave function corrections are included. 

Dynamical processes in a first order phase transition have been studied before with lattice 



simulations in the 3-dimensional Ising model ||25|. However, in these studies the parameters 



of the simulations were chosen so that the nucleation rate was relatively large, so that the 
nucleation timescale was, at most, only few orders of magnitude larger than the microscopic 
interaction timescale. Thus, the nucleation process could be observed simply by waiting for 
the nucleation from a metastable to a stable state to happen during a straightforward stan- 
dard (real-time) simulation. This is also the case for the recent work of Borsanyi et. al. [EB 



However, in this work we are interested in extremely strongly suppressed nucleation. 
Indeed, very slow nucleation rates are a quite common characteristic of the first order phase 
transitions in Nature: fundamentally, this is due to the fact that the external parameters 
which drive the transition (temperature, say) usually vary on several orders of magnitude 
longer timescales than the microscopic interaction scale. This physical situation is out of 
reach of any numerical real-time simulation method relying on spontaneous appearance of 
bubbles in the metastable phase. On the other hand, the method presented in this work can 
be applied to an almost arbitrarily slow nucleation. This method has also been used with 
the 3-dimensional cubic anisotropy model p7 . 



An outline of the paper is as follows. In Section |I|, we present the approach and discuss 
the obstacles. The discussion does not rely in any way on the specifics of the electroweak 



nucleation problem, except that it can be considered as a problem in classical statistical 
mechanics (both in the thermodynamics and in the dynamics of the system) . Subsections 
ll7\| and [11 B| , together with subsection [IV B| , are the most important parts of the paper 



to understand; we encourage the reader to concentrate on them. In Section |T|, we review 
why it is true that the physics of the electroweak phase transition can be considered, both 
thermodynamically and dynamically, as a classical statistical mechanics problem. This 
section is a review of previous literature, included mostly to make the paper self-contained. 



In Section |V| we present the numerical tools we use and the details of the calculation. It 
ends with a presentation of our results. Section |V| presents a number of alternative, "more 
perturbative" and less numerically intensive ways to try to determine the bubble nucleation 
rate, most of which have previously been used in the literature. We systematically compare 
these approaches to the nonperturbatively determined result, to analyze the reliability of 
the other approaches. Of these, the most trustworthy are the thin wall approximation using 
nonperturbative inputs (surface tension and latent heat), and two loop perturbation theory 
including Higgs field wave function corrections. Each of these approaches makes errors of 
order 20%; other approaches, including those most widely used in the literature, give results 



off by almost a factor of 2. Finally, Section VI presents our conclusions 



II. STRATEGY TO DETERMINE NUCLEATION RATE 

The next section will show why the calculation of the bubble nucleation rate at the 
electroweak phase transition is a problem in classical statistical mechanics. In this section 
we will assume this to be the case, and discuss the strategy for solving this statistical 
mechanics problem. The basic idea is that nucleation from the metastable symmetric phase 
to the stable broken phase is limited by the rarity of "critical bubble" configurations which 
lie in between. A "critical bubble" is, roughly, a configuration which in the medium termQ is 
about equally likely to evolve towards the broken and homogeneous symmetric phases. An 
outline of the strategy is: 

1. choose a measurable which will distinguish which field configurations are near the 
metastable symmetric phase, which are near the broken phase, and which are near the 
critical bubble; 

2. evaluate the probability to be in the (exponentially rare) critical bubble configurations; 

3. determine how quickly a "critical bubble" evolves towards one of the (meta) stable 
phases; 



We will distinguish three time scales. The short time scale is the longest time scale of typical 
thermalization processes in either phase, ~ l/g'^T. The long time scale is the time scale for 
nucleation to occur, ~ e^^^ /T. By the medium term we mean on a time scale well separated from 
the short and long time scales. At many points in our discussion there will be ambiguities of order 
the ratio of two of these time scales; but for strongly exponentially suppressed nucleation problems 
such ambiguities are tiny. 



4. determine the "dynamical pref actor," which tells what fraction of imputed "critical 
bubbles" really represent midpoints on a trajectory carrying the metastable symmetric 
phase to the stable broken one. 

This section elucidates what we mean by each of the above, and why the whole approach 



is possible. Our strategy is similar to hanger's classic method [^, except that the saddle 
point treatment of the "critical" configurations is replaced by a Monte Carlo calculation, 
and we take care to treat completely the microscopic dynamics during "barrier crossing." 
Essentially the same technique has already been applied to determine the broken phase 
"sphaleron rate" in [^9[^, but we will not assume the reader is familiar with those papers. 

Everything we say in this section is generic to first order phase transitions of liquid-gas 
type, ie. where there is no breaking of a global symmetry but the phases can be distinguished 
by the value of the volume average of a single, scalar local measurable. (However, this 
method is also fully applicable to transitions exhibiting a real global symmetry breaking.) 
In our case the volume averaged observable will be the average of the (length) of the Higgs 
field, (pl^ = (1/^) / c?^x2<l>^$. If it makes the reader more comfortable she can think of water 
and gas below the critical temperature, with the density as a measurable distinguishing the 
phases and pressure as a control parameter, or of a ferromagnet below the Curie temperature, 
with the magnetization as a measurable and an applied external field as a control parameter, 
rather than the electroweak model with 0^^ as a measurable and temperature as a control 
parameter. The big difference from the liquid gas system is that we know how to simulate 
the microscopic physics accurately even away from the second order endpoint, so Monte 
Carlo techniques applied to a first principles microscopic description of the thermodynamics 
can give reliable results. 

With suitable generalizations the method described here can be applied to almost any 
metastable state decay problem, provided that 

1. both the thermodynamics and the real time evolution of the system are amenable to 
numerical analysis, 

2. there is an observable which can unambiguously distinguish the phases and can resolve 
the potential barrier between the them, and 

3. the potential barrier is large and the tunnelling rate is small. 

The first two conditions are quite generic, and if the third condition is not satisfied (for exam- 
ple, in "quenching" type problems), then one can make straightforward real-time simulation 
of the decay without relying on the special methods described here. 

A. General picture of homogeneous nucleation after weak supercooling 

Fix the ratio of the scalar self-coupling and gauge couplings, X/g"^, to a value where there 
is a first order electroweak phase transition (or fix T to a value below the Curie temperature 
for a ferromagnet, or T below the critical temperature for a liquid-gas system), and ask 
how the canonical ensemble is distributed close to the critical temperature. In particular, 
consider how the constrained free energy F{(j)l^) varies as a function of 
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FIG. 1. Cartoon of how the constrained free energy = —log (probabihty of (/)^v) varies with 
at the equilibrium temperature in a large volume. The vertical axis gives minus the log of the 
fraction of states in the canonical ensemble with the given value of cj)^^. The dotted line gives the 
free energy of a spatially homogeneous configuration with that value of 0^^; the truly most probable 
configurations at intermediate values are mixed phase configurations. 

where the counterterm is needed to subtract UV divergences so the whole is well defined. 
We could consider any other measurable which has UV finite variance, but 0^^ will prove 
particularly convenient below. In a very large volume, at Teq (or zero external magnetic field 
in a ferromagnet, or the boiling pressure in a liquid-gas system), the constrained free energy 
density dependence on 0^^ will qualitatively resemble that of Fig. [^. The constrained free 
energy means —T times the log of the weight of configurations in the canonical ensemble 
which have the specific value of 0^^. Very low or very high values of 0^^ are extremely rare; 
but values intermediate between the two bulk phases have free energy per volume equal 
to those of the bulk phases. This is because, besides the "expensive way" of getting an 
intermediate value of 0^^ — having $^$ equal the desired value homogeneously through the 
volume of interest — there is a "cheaper way," which is to have part of the volume be in one 
phase and the rest in the other phase. Then the disfavored intermediate values of $^$ are 
only achieved in an interface between the regions, whose volume does not scale extensively 
with the system volume. (Note that the fact that intermediate values of $^$ are disfavored, 
is exactly the statement that there is a first order phase transition with different values of 
$^$ in the two phases. )Q 

What if we ask about a smaller volume, where the amount of space in the interface 



^Strictly speaking, when we talk about spatial variations of ^^^ or spatially homogeneous <1>^$ 
we actually mean a "coarse-grained" quantity: normally ^^^ fiuctuates wildly from point to point 
even in the pure symmetric or broken phases (it is, indeed, UV divergent). The coarse-grained 
$'1^$ is obtained by averaging over length scales ~ bulk correlation length ^. This is equivalent to 
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FIG. 2. Same as Fig. |I| but in a finite volume where the free energy of the interface between 
phases is not considered neghgible. The free energy of a mixed phase state is higher than either 
pure phase because of the surface tension of the phase boundary. The figure also illustrates the 
physical appearances of the states which dominate the ensemble at intermediate values of (/>^y. 

between extensive phases is not negligible? A qualitative cartoon of the answer is given in 
Fig. ^. There is a free energy "barrier" between the two phases due to the free energy cost 
of the interface separating the two phases; but it is much lower than it would be if we had 
to stick with spatially homogeneous intermediate states. The barrier is roughly the surface 
tension of the interface times its area, which scales as the length squared of the box, while 
an extensive quantity would scale as length cubed. 

The free energy at values of 0^^ between the stable phases gives us information about 
the free energy cost of mixed phase configurations; in particular, the free energy near the 
symmetric phase tells about the cost to have a small bubble of the broken phase in the 
symmetric phase. To see the value of this, in determining the rate of bubble nucleation, we 
now discuss how the picture changes when we change the temperature. The short answer 
is that one should "tip" Fig. ||, adding a linear in <^\^ term to the free energy. In fact, for a 
special measurable this statement is exact, as we now discuss in some detail. 

In the 3-D effective theory approximation we are working in (see next section), a variation 
bT of the temperature corresponds to a change Sm'^^rp in the thermal Higgs mass squared. 



The size of the change can be read off from Eq. ( p.3| ). Henceforth we will only talk about 
changing m^jj,. The way that one determines the constrained free energy plots we have been 
discussing is that one finds the probability that a configuration drawn from the canonical 
ensemble has the value of 0^^ of interest. If we are interested in the free energy as a function 



integrating out spatial momenta larger than ^ ^. In the pure bulk phases the coarse-grained $^$ 
is almost homogeneous by construction, and non-homogeneous only in the mixed phase. 



of an operator O, we want to know 

li^ = -\n J ViA,<!>)expi-H/T) 5(o(A<|.)-Oo), (2.2) 

up to an overall constant which is uninteresting (and depends on how we define the normal- 
ization of the path integral). For the special case that our operator is 0^^, the constrained 
free energy has an extremely convenient property. Observe from Eq. ( p.7|) that the way that 
m'j^rp enters the Hamiltonian is 

H = Hm=mo + I^il^ J d'x2<^^^x) , (2.3) 
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but we can now use the delta function to replace 0^^ in the exponential with 4>1^q, which is 
not integrated over; pulling it out of the integral gives 
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The second term here is independent of m^uT'i i^ can be determined once and used at any 
value of m^HT thereafter. Hence the effect on -F(0av) of shifting rn^rp is very simple; it just 
adds an extensive, linear in 0^^ term to F. (The same thing would happen if we considered 
the magnetization in a ferromagnetic system where we vary the external magnetic field, or 
the density in a liquid-gas system where we vary the pressure. The key is to consider a 
measurable which appears in the Hamiltonian next to the control parameter which takes us 
through the transition.) If we used a different measurable the qualitative behavior would 
be the same-the free energy as a function of that measurable would be roughly a "tilted" 
version of its Teq appearance-but this would not hold as an exact quantitative statement. 

Now consider how the free energy plot looks when we shift m\jrp. A cartoon is provided 
in Fig. ^ For small amounts of supercooling, the symmetric phase minimum shifts over 
slightly, but persists as a local minimum of the free energy. It is labeled {A) in the figure. 
Since we will be concerned with classical dynamics, for an element of the thermal ensemble 
at {A) to get to the global minimum, it must pass through (5), (C), and {D). The time 
evolution of a configuration at {B) is almost certain, in the medium term, not to go to (C), 
because F/T is — log(Probability); there are vastly more states with 0^^ equal the value at 
[B) than the value at (C), so only a tiny fraction will evolve to (C) in the medium term, 
since time evolution preserves the canonical ensemble. On the other hand, configurations at 
{D) are almost certain not to "go back," and will continue to the broken phase minimum. 

8 



Free energy 



Free Energy Plot at Equilibrium 
Temperature 



Free energy 



Magnified Region 



Free Energy Plot Below Equilibrium 
Temperature 




-> 




A: metastable minimum 
B: "too small" bubble 
C: "critical" bubble 
D: "too big" bubble 



FIG. 3. Cartoon showing how free energy in a finite box changes when we lower the tem- 
perature. For small temperature changes, both minima survive, but one (A) is no longer globally 
stable. The least likely configuration on the way to the stable minimum is (C) the critical bubble: 
its unlikelihood restrains the rate at which configurations near (A) go to the true minimum. 

It is the rarity of configurations at (C) which limits the rate at which configurations near 
(A) time evolve into the broken phase. If it were true that every configuration at (C) were 
on the way from {A) to the broken phase (or going from the broken phase to (A)), then 
with a little dynamical information we could determine the nucleation rate. It may be, 
though, that most configurations at (C) are either coming from (B) and going back there, 
or coming from {D) and going back there. (This might happen if the choice of measurable 
is not optimal, for instance.) In this case, nucleations from the symmetric to the broken 
phase would be even rarer than the free energy of states at (C) implies. So the rarity of 
configurations at (C) provides an upper bound on the nucleation rate, which can be turned 
into a determination with some additional dynamical information. 

B. Real time rate, dynamical prefactor 



Now we discuss how to turn the discussion and cartoons of the last section into a calcula- 
tion of a real time rate for nucleations. Suppose we have, by multicanonical tools discussed 
in subsection [IV B| , computed the constrained free energy as a function of some measurable, 
which we take to be 0^^. It is also possible for us to collect a sample of the canonical en- 
semble restricted to some narrow range of (pl^, for instance, the range right around (/'avC! 
the least likely value of 0^^. What to we do with them? 



The first thing to do is to determine how exponentially suppressed critical bubbles are. 
But the answer depends on the measurable and does not give a real time rate. The second 
step is to determine what the flux of states in the canonical ensemble through the critical 
bubble is, normalized to the probability to be in the symmetric phase. That is, we should 
determine 
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with P (condition) denoting the fraction of the canonical ensemble which satisfles the con- 
dition, and with e inflnitesimal. The flrst term here is the probability density to be at the 



critical value 
be written as 



^L,c of 



which is our "deflnition" of the critical bubble. It could equally 
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and will be evaluated by Monte Carlo. The second term in Eq. ( |2.6| ) is the mean change, in 
absolute magnitude, of (pl^ in a time interval At, sampled over configurations at the critical 
bubble. Provided we take At shorter than any typical infrared scale — in particular it 
should be shorter than the time scale to go from a configuration with 0^^ = 0^^ ^ to one 
where F/T differs by order 1 — then the ratio will tell the number of configurations which 
pass from one side of the critical bubble to the other in a time interval At, times the time 
interval and divided by the number of symmetric phase configurations. In other words, the 
combination gives the fiux of configurations in the canonical ensemble through the critical 
bubble. 

This fiux is A^OTthe bubble nucleation rate we are after, though it is clearly an upper 
bound. We have to multiply by the fraction of critical bubble crossings which actually 
mediate a change from phase to phase. To this end we define a "dynamical prefactor," d. 



as 



d = 



trajectories getting from (B) to (D) 
crossings of (C) 



(2.8) 



In other words, if we consider all real time trajectories, d is the fraction of crossings of (C) 
which represent "permanent changes" from one side to the other of the barrier. To determine 
it we sample the ensemble of configurations restricted to those at the critical bubble, and 
for each element of this ensemble we construct a trajectory forward and backwards in time, 
long enough to see the configuration come from and go to an exponentially more common 
value of 0^ . Sampling trajectories every At, d is 



1 



# of crossings 



X 



change sides 
don't 



(2.9) 



The average is over the canonical ensemble, restricted to configurations with 10^^ — 0avcl < 
e/2, and over trajectories through those configurations. The measure to be used is the 
canonical one, times |A0^y/At| evaluated where the crossing takes place; so the sample is 
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precisely the sample of the flux of states through 0^^ = 0^^ q. In Eq. (pl9|), (# of crossings) 
means the number of crossings of the critical bubble a trajectory makes in the medium term. 
To determine this we need to follow the trajectory until it reaches an exponentially more 
common value of 0^^, ie. point (B) or (D) in Fig. ^J^ The trajectory "changes sides" if 
it goes from {B) to (D) or vice versa, and does not change sides if it returns to the same 
side it came from. In our case, the real time evolution will be Langevin evolution, see next 
section, and the forward and backwards time evolutions are just two Langevin evolutions 
with two different realizations of the random force. For the case of Hamiltonian dynamics, 
the forward evolution would be evolution with a set of momenta drawn from the canonical 
ensemble and the backwards evolution would be evolution with the same momenta, but sign 
reversed (T conjugated). In either case we are approximating the expectation values in Eq. 
( p.9| ) with an average over a sample of trajectories, ie. we take the average in Eq. ( pT9| ) by 
a Monte Carlo integration. We discuss why this procedure gives the correct nucleation rate 
at more length in Appendix ^. 

For Hamiltonian dynamics, in an UV regulated theory, both d and (|A0^^/At|) should 
have well deflned small At limits. This is not the case for Langevin dynamics, however. If 
we sample a Langevin trajectory with a smaller At, the number of crossing should grow, as 
each crossing gets "resolved" into potentially more; this is a normal feature of a Brownian 
path. However, (|A0^y/At|) will also depend on At by a compensating amount. For suitably 
short At, the time history of 0^^ near each crossing looks like a Brownian random walk. By 
well known properties of Brownian random walks, (|A0^y/At|) scales as (At)^^^"^, while d 
scales as (At)^/^, and the product has a flnite small At limit. Hence, for Langevin dynamics, 
neither d nor (|A0^^/At|) are well deflned but the product is. 

It is flnally the product, 

nucleation rate 1 , , .,. r, , /^..^n 

= — probability flux x d, (2.10) 

Volume 2V 

which we are interested in. The factor of (1/2) is because half of the permanent crossings 
the algorithm flnds are into the symmetric phase. The factor 1/V turns the nucleation rate 
into a rate per unit volume. 

C. Complications: peculiar behavior in finite volumes 

In this subsection we discuss some complications with applying our technique, arising 
from flnite volume effects. The conclusion will be that the volume must be fairly large, so 
the bubble interior flUs at most about 15% of it; and that as a consequence, it is best to 
choose a measurable with a very small variance in the metastable phase. The impatient 
reader may want to skip this section and just accept that conclusion. 



^We must also check that this criterion is sufficient to ensure that it is very unlikely for the 
trajectory to return again to (C), which is not necessarily ensured; the ensemble of configurations 
at {B) which have just evolved from configurations at (C) is not the same as the ensemble of all 
configurations at {B). In practice this does not prove to be a problem. 
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What we want is the nucleation rate in very large volumes. In practice it is impossible 
to work directly in a very large volume, for reasons of numerical cost. Naively the cost of 
performing the Monte Carlo calculation scales as the volume, but in practice the scaling 
is still more severe, because of memory and communication costs and because the Monte 
Carlo becomes less efficient in large volumes, particularly at the value of 0^^ where the typical 
configuration changes from being homogeneous, supercooled symmetric phase to being an 
isolated bubble. 

a. Resolving the critical bubble: There is another reason why we would like to work in a 
small volume. This is because any volume averaged measurable, for example 0^^, fiuctuates 
also in the pure symmetric or broken phase. The width of these fluctuations behave as 
1/vK, V the volume, whereas the contribution of a fixed size bubble to 0^^ scales as l/V . 
Thus, if we keep the size of the critical bubble constant but increase V , the fiuctuations 
of 0^v i^ the pure symmetric phase outside of the bubble increase and degrade the the 
"cleanliness" of 0^^ as a description of the critical bubble. 

Concretely, what we mean by this degradation is that the more symmetric phase there 
is, the more likely it is that the value of 0^^ consistent with the critical bubble (point C in 
Fig. 1^) really arises from too small a bubble plus an upward fiuctuation in the symmetric 
phase contribution to 0^^, or too big a bubble and a downward fiuctuation. The result is that, 
as we increase the volume, the measured free energy of the bubble should fall somewhat faster 
than the "true" bubble free energy (which decreases as — log(l^), due to the translational 
zero modes of the bubble configuration). However, the measured value of the "dynamical 
prefactor" d, Eq. ( p. 91) , should become smaller, since it is more likely that an imputed 
critical bubble is really on one or the other side and will begin and end at the same phase. 
These effects cancel exactly, so that the nucleation rate has a well-defined infinite volume 
limit. However, as d becomes smaller, it takes more work to measure it with good relative 
accuracy. This effect also degrades the efficiency of the Monte Carlo. 

b. Maximum size of the critical bubble: For the above reasons it is to our advantage to 
use as small a volume for a given size bubble as we can get away with. Naively, this means 
we should use the smallest volume for which the critical bubble is unable to "see itself" 
around our periodic boundary conditions (we choose to work in a cubic box with periodic 
boundary conditions. It appears necessary to use a box which has everywhere a fiat spatial 
metric, and does not have boundaries, because either effect could modify the bubble free 
energy.). Very naively, this means we must ensure that the radius of the critical bubble rmax 
is less than half the box length L, r^^^ < L/2. 

In fact the criterion for a sufficiently large L is more severe than L > 2rinax- The reason 
is that a spherical broken phase bubble with L/3 < r < L/2 is at best only metastable: a 
configuration with a cylinder of the broken phase, extending through the length of the box 
and having the same volume as the bubble, will have smaller phase interface area and hence 
smaller total free energy. To compare the favorability of different geometries as a function 
of broken phase volume fraction, we will make a "thin wall" approximation in which the 
interface between phases is treated as a geometrical surface, and the free energy is equal 
to its area times the surface tension: F = a A. This approximation is correct in the limit 
that the box size L is much larger than the wall thickness. Though our simulations will 
not be strictly in this limit, it is a suitable approximation for understanding the relative 
favorabilities of different mixed phase geometries. One can then write down how the area 
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Geometry 


Sphere 


Cylinder 


Planes 


Area(r) 


47rr2 


2'KrL 


2L^ 


Volume (r) 


47r 3 
3 ' 


Txr^L 


any 


A{V) 


(367ry2)i/3 


2^/txLV 


2L2 


Volumes where stable 


- 47r/81 


47r/81 - I/tt 


lA - (i-iA) 


Vol. where metastable 


- 7r/6 


l/47r - 7r/4 


- 1 



TABLE I. Area and volume as functions of r, and area as function of volume, for each possible 
geometry for phase coexistence; and derived volume range where the geometry is preferred and 
where it is metastable in the strict thin wall limit. The upper end of the sphere and cylinder 
metastability ranges are where they touch themselves across the periodic boundary conditions. 
The lower metastability limit on the cylinder is where it becomes unstable to a (5r = sin(27rz/L) 
excitation. 

and volume vary as a function of radius for a sphere and a cylinder extending the length of 
the box, and for a pair of planar interfaces. The results are shown in Table | and compared 
in Figure §. Equating the areas at fixed volume. 



Aiy) 



sphere 



(SevrV 



2N1/3 



V47rLr = A{y) 



eylinder ; 



(2.11) 



we find that the sphere and cylinder are degenerate for V^ = L^ x 471/81 ~ 0.155L^, and the 
cylinder and planes are degenerate for V = L^/tt. A plot of the free energy in a very large 
box should approximately resemble the solid line in Fig. §. It will have slope discontinuities 
where one geometry gives way to the next. The metastability of one geometry against 
another can be quite strong, and this makes it challenging to perform a multicanonical 
Monte Carlo determination of the free energy in a volume large enough that the interfaces 
are thin compared to L, unless one only wants the plot in the region where one geometry is 
relevant. (Fortunately this will be the case for us.) 

Since both the cylinder and the planar geometry have free energies which clearly depend 
essentially on the box geometry and volume, the part of the free energy plot where they 
dominate describes finite volume artifacts rather than physics which has any correspondence 
to that at larger volumes. However, the radius of the sphere where the cylinder becomes 
equally favorable is only r = L/3, safely small enough that the sphere will not "see itself" 
around the periodic boundary conditions (unless L is only a few times the interface width). 
Hence we can expect the sphere geometry, until it stops being the favored geometry, not 
to care about the periodic boundary conditions but to represent more or less faithfully the 
behavior of an isolated bubble in a larger volume. 

We may also expect that, in a suitably large volume, the transition rate between the 
spherical and cylindrical geometries is sufficiently slow that we could rely on the metastability 
of the sphere to explore larger radii than r = L/3. However we will conservatively not do 
so in what follows, but will restrict ourselves to such volume and bubble size combinations 
that the critical bubble we obtain will have 0^^ < 0.150^^ (broken) + O.850^y(symmetric). 
We will also check to see that the bubbles we analyze are approximately spherical and not 
cyUndrical and that they do not touch across the periodic boundary. 
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0.4 0.6 

Volume fraction 

FIG. 4. Area, and hence free energy, as a function of volume fraction, in the thin wall 
approximation. The solid line is the minimum over interface geometries; the large volume free 
energy curve would follow the solid line. Dotted lines are the metastable extensions of the sphere 
geometry (sloping) or the planar boundary geometry (flat), while the dashed lines show metastable 
extensions of the cylindrical geometry. 

III. ELECTROWEAK BUBBLE NUCLEATION AS A PROBLEM IN CLASSICAL 

STATISTICAL MECHANICS 

In this section we briefly review why we can view the bubble nucleation problem as a 
problem in classical statistical mechanics for Yang-Mills Higgs field theory. Nothing in this 
section is new; it reviews the last few years' developments both in the thermodynamics of 
the electroweak phase transition and in the dynamics of infrared Yang-Mills Higgs fields. 
We include it here to make the paper more self-contained. Readers who are already familiar 
with this material may want to skip this section. 

A. Thermodynamics: dimensional reduction 

Here we review how the thermodynamics of infrared fields in the SU(2) sector of the 
standard model is well approximated by a 3 dimensional path integral, which is the same as 
the partition function of classical 3+1 dimensional SU(2) Higgs theory at finite temperature 
(and with suitable regulation and counterterms) . 

The full theory we are interested in has thermodynamics described by the path integral 
(Lorentz indices are Euclidean with positive metric; Greek indices range over space and 
time, Latin indices i,j, k only over space, indices a, b, c are SU(2) group indices) 
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^1 



Z = J P(<l>, A^, etc) exp{-SE/h) 

1 



h/T 



dr / d'^x 



—.F^^^Fr + {D,<l>y{D^<l>) + (mL + A<l>t$)$t$ 



(3.1) 



+ hypercharge + fermions + glue 



Here and throughout g is g^, the weak couphng. The r integration has periodic boundary 
conditions for bosons and antiperiodic boundary conditions for fermions. Beginning here 



we will neglect hypercharge, to simplify things. This is not too bad an approximation ^1 
Including it would be a straightforward extension of what we discuss. 

There are two things to observe right away about this theory. First, mean field theory 
predicts that if mjj^ changes, there is a second order phase transition. Second, at the scale 
T (which will be of order the weak scale, T ~ 80GeV), the coupling is weak. This is just 
the statement that the weak sector of the standard model is indeed weakly coupled. Hence, 
if there is a phase transition, barring some large hierarchy of couplings such as X/g'^ <^ 1, it 
will be weak, and correlation lengths will be ^ ^ 1/T. 




thickness=l/T 



FIG. 5. Cartoon of how a 3+1 dimensional spacetime, K^ x S"^, drawn here as 2+1 dimensional, 
can look effectively 3 dimensional for long distances. 

As motivated in Fig. ^, at such infrared scales the effective behavior is 3 dimensional. 
This is just because any field varying only on the length scale ^ 3> 1/T will not vary ap- 
preciably across the Euclidean time width of the "slab." One sees this formally by Fourier 
transforming the r direction in Eq. (|3.1| ). The Euclidean frequencies arise from transform- 
ing a compact range and so are discrete: (9t-$(t))^ becomes (27mT)^$^ for bosons, while 
?/'7°9t-^(t) becomes ((2n + l)7rT)?/'„7o^„ for fermions. All but the n = bosonic mode are 
very heavy and can be integrated out. The result (continuing to use 4 dimensional notation 
for fields and couplings) is 



Se 

h 



1 

T 



d'x ^Ft^F^^ + ^{D^AoTiD^Aor + (A$)^(A$) + m^r$t$ + A($t$)2 
+m^ASAg + i^_t^^i^^A«^^$t$ + 0((72a^)(A^A^)2 + Dim 6 



(3.2) 



where "Dim 6" indicates dimension 6 and higher induced operators, which have an irrel- 
evantly small effect on the infrared physics, so we immediately drop them. The fields we 
write here correspond to the zero frequency components of the 4 dimensional fields, ie. $ in 
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Eq. ( p.2|) is T J dT^{x, r), with $(x, r) the field appearing in the path integral in Eq. (p.ll). 
In this expression A and g^ are couplings of an effective 3D theory, which (after dropping 
dimension 6 operators) is super-renormalizable; they do not run. Their relation to the coef- 
ficients of the full theory, and a detailed discussion of the matching procedure used to derive 
them, is given in |0. In particular we mention that the Higgs mass squared, rn^jj,^ receives 
O^g^T"^) positive thermal corrections; 

mlr = mL,en + ^^!±M±^T^ + ( /. dependent 0{g'T^)) . (3.3) 

Here yt is the top quark Yukawa coupling. Since mvac,ren is negative, varying T can change 
the sign of the Higgs mass squared and induce a phase transition.]] For generic values of 
X/g"^, correlation lengths at the phase transition are of order C, ~ 1/g'^T, which is why we 
can drop the dimension 6 operators. 

Note that the form of S^/h looks very much like the Boltzmann factor, E/T, of a classical 
theory. The main difference is that in classical SU(2) Higgs theory, we expect the electric 
field strength and the Higgs field momentum H, rather than Ao, to appear; 



H = I d^x 



-^F^^F'J + \EtEt + H^H + (A$)^(A$) + m]jT^^^ + \{^^^f 
Ag^^ ■' 1 



(3.4) 



Here the delta function enforces Gauss' Law. But as noted by Ambj0rn and Krasnitz ||32|| , 
if we implement Gauss' Law by introducing a Lagrange multiplier, which we suggestively 
name Aq, 



bi... 



J VAo exp [i J d'xA^, {{DiEif + [$^^n + c.c] ) /t) , (3.5) 



then the E and H integrations are Gaussian and can be performed, generating 

1 0^ 

H D -{D^AoY{D^AqY + ^AlAl^^^ + (0 + radiatively induced )AIAI . 



(3.6) 



Here the bare Debye mass is zero but one is radiatively induced, with a linear divergent 
coefficient in any regulation where such terms do not identically vanish (such as the lattice). 
Hence, the thermodynamics of the full theory DO look like those of the classical theory, 
except that the Debye mass of the classical theory is radiatively induced and regulation 
dependent, and there are very small 0{g'^aw) extra interaction terms involving the Aq field, 
present in the actual thermodynamics but not in the thermodynamics of the classical system. 
It is a good approximation, for the thermodynamics both of the classical theory and of the 
dimensionally reduced full theory, to integrate out the Aq field, including it by the radiative 



^Note that both m\j< and mj^ ~ g^T'^ renormalize logarithmically at the two loop level, 
{dm'jjrp / d In fi) ~ g'^T'^, so the sign of m^i^rp near the transition is actually renormalization point 
dependent. 
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corrections it will induce in the remaining couplings. In this approximation one shifts slightly 
the coefficients of the terms in Eq. (|3.2|) which do not contain Aq, and drop those which do. 
The change in the coefficients is computed in 0. After this approximation, the classical 
thermodynamics and the dimensionally reduced thermodynamics coincide exactly. The final 
form of the partition function describing the thermodynamics is then 



Z = Jv{A,<!>)exp{-H/T), 

(3.7) 



H = d'x 



1 
4^ 



ir^Wtj + (A$)^(A$) + m^y$t$ + \{^^^f 



If this final integration over the Aq field is not deemed reliable enough it is straightforward 
to modify what we will do below to include it in the thermodynamic calculation. 

We will also mention a few results, perturbative and nonperturbative, which have been 
obtained for the partition function shown above. Perturbatively we can describe the strength 
of the phase transition by studying the effective potential for the (gauge fixed) Higgs field 



'2<l>jQ^^$cond- Note that there is no good gauge invariant, nonperturbative way 
to separate the condensate from the fluctuations. When gauge field fluctuations become 
large, the whole perturbative approach becomes questionable. Nevertheless, when the phase 
transition is strong, perturbation theory is useful, essentially because gauge field fluctuations 
are suppressed in the broken phase (which is therefore well described). At one loop, and 
neglecting scalar loops as is appropriate in the Xj g"^ <^ 1 approximation, 

^iioop = ^0^-^0^T+^0^ (3.8) 

Z iOTT 4 

This effective potential can have two minima because of the (one loop) negative (f)^ term. 
Since it is a loop effect which allows a ffist order phase transition, we say the transition is 
radiatively induced ffist order. At the transition, the broken phase value of is such that 
the cubic term and (tree level) Xcj)^ term are of the same order. Since this requires a one 
loop effect to be of the same size as a tree level one, it either implies that \/ g'^ <C 1, or 
that perturbation theory will not be a reliable expansion. Hence perturbation theory deter- 
mines attributes of the transition at best as an expansion in \/ g'^. It is this relatively poor 
performance for perturbation theory which makes a nonperturbative treatment necessary. 

Nonperturbatively it is known that, as expected, the phase transition is well described 
by perturbation theory for small \/ g"^'-, but perturbation theory is completely wrong for 
larger values 0J31- ■'■^ ^^^^ there is NO phase transition in the MSM above a critical value 
\/ g'^ ~ 0.0983 ± 0.0015 ^j. In extensions of the standard model with new bosons which are 
light at the phase transition, we must include the new light bosons in the effective theory 
considered. At least for the case of an added scalar top, the strength of the phase transition is 
significantly enhanced [p3| , p^ . It would be straightforward but more numerically expensive 
to apply the tools developed here to this physically interesting case. 

B. Dynamics: classical effective theories 

The infrared thermodynamics of the SU(2) sector of the standard model match those of 
classical 3+1 dimensional SU(2) Higgs theory, as discussed above. Does this matching also 
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apply at the dynamical level, as originally conjectured (in a 1+1 dimensional context) by 
Grigoriev and Rubakov [^|? In other words, are the dynamics of the infrared SU(2) Higgs 
fields described by classical Hamiltonian dynamics? 

The answer is "no, it is more complicated than that."0 As we now discuss, the dynamics 
of infrared gauge and Higgs fields are indeed classical; but they are not described by classical 
Hamiltonian dynamics. They are, to leading order in the logarithm of g, described by 
"classical" Langevin dynamics. Because a factor of 2 error in the treatment of the dynamics 
will change our determined rate by ~ exp(±l), while the rate itself is ~ exp(— 100)T^, we 
will find it sufficient to make this approximation and take the dynamics to be Langevin. 

To see that the dynamics of the real system are not classical Hamiltonian dynamics, look 
again at the thermodynamic discussion, especially Eq. ( ^^ and Eq. ( p.6|) . We see that there 
are linear UV divergences in the classical Debye mass for the Aq field; and yet the Ao field 
is related to the electric fields E, which generate the Hamiltonian dynamics for the gauge 
fields. This implies that there are divergent UV corrections to the gauge field dynamics of 
the classical fields. As first advocated by Bodeker, McLerran, and Smilga ||36|, we should 



consider this a potential problem for the study of the classical field dynamics. In fact, as 
argued by Arnold, Son, and Yaffe [^, what it means is that the classical Hamiltonian 



dynamics do not have a good regulation independent limit. Hence, the infrared gauge field 
dynamics of the classical theory technically do not exist. The dynamics of the full theory 
can scarcely coincide with those of the classical theory if the classical theory's dynamics are 
sick. 

The physical origin of this problem is actually well known plasma physics. Transverse 
electric fields in a plasma feel Landau damping. This leads to very slow, overdamped 
evolution of infrared magnetic fields, as is typical in a conducting medium. As the classical 
theory cutoff is lifted, there are more and more "plasma" degrees of freedom, and the 
damping becomes ever more efficient. The correct treatment is to make the damping have 
the same efficiency as in the quantum theory. This requires studying the classical theory with 



hard thermal loop (HTL) effects I^H] included. (The hard thermal loops are the nonabelian 
generalization of Debye screening. Landau damping, and other plasma effects familiar from 
electromagnetic plasmas.) Such a classical, but HTL including, treatment should be correct 
at leading order in the coupling g. Two numerical implementations of such a classical theory 
now exist; one |^ is based on a proposal by Hu and Miiller |^, and one [|17| is based on 



a proposal by Bodeker, McLerran, and Smilga |3^, and more recently discussed by lancu 
PT| . Both are extremely complicated. Probably the second method could be utilized in the 
type of computation we are going to discuss, but the numerical effort would be substantially 
greater than what we discuss below. 



As first demonstrated by Bodeker |[T5|, and further discussed and clarified both by 



Bodeker ^^, Arnold, Son, and Yaffe [§^, and Litim and Manuel |0], inclusion of the 



In defense of Grigoriev and Rubakov we should mention that the complications discussed in 
this subsection do not arise in 1+1 dimensions, where the UV behavior is much more mild; hence 
their conjecture is correct for the problem they were addressing, namely the dynamics of the 1+1 
dimensional abelian Higgs model. The problem is applying it to the 3+1 dimensional problem of 
interest instead. 



hard thermal loops in the infrared equations of motion is actually unnecessary, because at 
leading order in log(l/5') the dynamics of the gauge fields are simple Langevin dynamics. 
We will not attempt to reproduce their arguments in detail here, but will only physically 
motivate them. 

First, consider the behavior of an abelian plasma. We will distinguish two characteristic 
length scales: the Debye length /d ~ I/'^d ~ ^/qT, and the scattering length /gcatt ~ 
1 / g^T log{l / g) . The former is the shortest length scale where plasma effects are important. 
The latter is the mean length over which a current can propagate, before it is disrupted by 
collisions in the plasma. It coincides with the free path of a charge carrier to undergo large 
angle scattering. 

On scales much longer than /scatt, a magnetic field evolves as if it were in a conductor: 
D X B = j = (Jc\E = a^iDoA, where o"ei is the electric conductivity. This is just Langevin 
dynamics; the time derivative of the field DqA is proportional to —dH/dA = D x B. In 
a thermal bath there is also a "noise" term uniquely determined by the thermodynamics 
(fluctuation dissipation). On scales between /^ and /scatt the story is significantly more 
complicated because the homogeneity scale is shorter than the free path over which an 
electric current propagates. In this regime the field evolution is described by an equation 
with a wave number dependent conductivity, which is nonlocal in real space. 

Now consider the case of a nonabelian theory. The key difference is that, for the non- 
abelian theory, /gcatt ~ ^/9'^Tlog{l/g). The reason is that nonabelian collisions exchange 
nonabelian charge ("color"), so any collision, however soft, can destroy the current a particle 
is carrying. Hence /scatt is of order the mean free path for any scattering, not just large angle 
scattering. Hence, on the length scale 1/g'^T, the dynamics are given by a simple Langevin 
equation 
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(3.9) 



with Nc = 2 for our SU(2) application. Here Uei is the nonabelian ("color") conductivity. 

The extension of these ideas to the case where there is a Higgs field turns out to be 
remarkably simple; one also evolves the Higgs fields under Langevin dynamics, but giving 
the Higgs fields a much larger diffusion constant than the gauge fields. For a discussion see 
These are the dynamics we will apply for the real time part of our studies below. Note 



that they are only justified at (next to |jT8|) leading order in l/log(l/f7), not a very good 
expansion; but we are willing to accept an approximation which will yield an error of ±1 in 
the exponent of the nucleation rate, since the rate itself is ~ exp(— 100). 



IV. COMPUTATIONAL DETAILS AND NUMERICS 

A. Our choice for a measurable 

In the discussion above we always chose to consider 0^^, the space averaged Higgs field 
length squared, as the observable used to distinguish the phases and the critical bubble. This 
has been the traditional "order parameter" observable in Monte Carlo simulations of SU(2) 
+ Higgs theory. It is easy to measure, and, because its variance is UV finite, given large 
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enough volume, it can unambiguously separate the symmetric and broken phases. Also, as 
we have seen, it is an extremely convenient choice because it makes it quite easy to use one 
set of multicanonical data to study a range of temperatures. Further, after UV counterterm 
subtractions, it has a good zero lattice spacing limit. 

However, a priori it is not obvious that this (or almost any other) measurable can distin- 
guish the critical bubble well enough for practical calculations. In particular we might worry 
that there is too large a "noise" contribution from the 85% of the volume which must be in 
the symmetric phase (see the discussion in subsection [II C|) . A necessary (but insufficient) 
criterion for a good measurable, needed to avoid this problem, is that it has a small variance 
in the symmetric phase. It turns out that 0^^ does indeed have a very small variance in 
the symmetric phase. The leading order perturbative result for the variance of 0^^ in the 
symmetric phase, in a volume V much larger than 1/m^ ^^, with rrisymm the symmetric 
phase scalar mass, is 
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(4.1) 

where in the last approximate equality we have substituted in the equilibrium, one loop 
symmetric phase Higgs mass for m. This variance is to be compared to the broken phase 
variance, which gets an added contribution from fluctuations in the zero mode, which has a 
condensate 00^ 

2 _ 202T 2 _ 4T 2 

(A?,, , broken ota "*"</>?,,, symm i\ /_9\_9t^ "*" (A?,,, symm ' \ I 



^2^^ ^"0L,symm [Xjgl^gly 



which is much larger for small [\j g^\ To get a phase transition strong enough to preserve 
baryon number after its conclusion, we will consider the case Xj g^ = 0.036; for this case 
the broken phase variance is about 100 times larger, and if 15% of the volume is in the 
broken phase, the variance contributed by this volume greatly exceeds that contributed by 
the symmetric phase; symmetric phase fluctuations will not pose a problem. Also note that 
the broken phase fluctuations are dominated by the motion of the condensate; and we expect 
that fluctuations in the condensate size are directly important to whether a bubble is more 
or less than critical, so these fluctuations may also not be dangerous. 

A posteriori we will of course determine whether the choice of measurable was a good one. 
For instance, when determining d, we can determine what fraction of trajectories crossing 
't'lv — '/'avc actually lead to a nucleation. We will flnd that 4>l^ is in fact a good measur- 
able; the fraction of trajectories crossing the critical bubble which lead to a nucleation is 
statistically compatible with 1/2, which is the maximum possible under Langevin dynamics. 
Note however that the above arguments suggest that, if we were studying nucleation out of 
the broken phase, then 0^^ would probably be a very bad measurable. We will not discuss 
what measurables might be useful for studying this nucleation rate, which fortunately is not 
cosmologically interesting. 

Note that the value for the scalar self-coupling we use, produces a Higgs mass lighter 
than the experimental limit; in fact there is no physical Higgs mass which gives such a small 
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value for the ratio of \/ g^ (parameters of a 3-D effective theory) ^ . We study this case as a 
toy example, because the phase transition is relatively strong here and perturbation theory 
arguably should be reasonable. 



B. Multicanonical method 

As described in Sect. |I|, by far the dominant factor in the bubble nucleation rate is 
determined by the constrained free energy F(0av) — ~^log-Pcan.(0av)' where -Pcan.(0av) i^ 
the canonical probability distribution of 0^^ for a particular value of rn^jrp-. 

Pcan.(0L)oc/P(A„<|.)exp[-i7/T] 5l^^l-Y.2^^^/v\ . (4.3) 

(The 2 next to $^$ is for the customary complex normalization of the Higgs field.) This 
probability distribution has to be determined for the whole range of values from the sym- 
metric phase to somewhat beyond the critical bubble value 0avC' see Fig. ^. 

In principle, -Pcan.(</'av) can be calculated with a standard lattice Monte Carlo computa- 
tion, where the configurations are sampled with the canonical probability 

Pcan. oc exp(-/f/T) , (4.4) 

where H is given in Eq. ( [3.71 ). Algorithms for such a sampling are well known, typically tak- 
ing the form of a Markov chain in which each configuration is a relatively small modification 
of the previous one. 

However, in the problem we are interested in, the probability can vary by a factor of 
~ exp(lOO) over the range of 0^^ of interest. A finite, canonical sample will simply contain 
no representatives for much of the range of 0^^ of interest, and hence give no information 
on the free energy in that part of the range; hence it will fail to determine Pcan.- For our 
problem, the canonical Monte Carlo method is utterly useless. This is exactly the kind of 
problem where multicanonical Monte Carlo methods excel. 

In a multicanonical simulation, the Monte Carlo sampling probability of configurations is 
modified so that the whole 0^^ range of interest is sampled with an approximately constant 
probability. This is achieved by sampling the configurations according to the probability 



Pmuca OC CXp 



H/T + W{<t>l) , (4.5) 



where the weight function W^(0av) is carefully tuned so that the multicanonical probability 
distribution 

Pmuca(0L) « / ^( A, ^) exp \-H/T + W] 5 Ul - Y: 2<l>t$/r^ 



oc exp 



wirj PcnXrj (4.6) 



is approximately constant. This condition is met if 



W{(j)l) ^ - logPean.(0av) + COUSt. (4.7) 
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The canonical expectation value of any observable O can then be obtained from a multi- 
canonically sampled set of configurations by reweighting the individual measurements with 
the weight function: 

(O) = Y, 0fce-^(<^-,.)/ ^ e-^(<^-,.) , (4.8) 

k k 

where the sums go over all configurations in the sample. It is not difficult to find an algo- 
rithm to perform the multicanonical update; if one has a Markovian canonical Monte Carlo 
algorithm, application of Metropolis accept-reject under the weight function exp(W {(pl-^)) 
after each update yields a multicanonical algorithm. As (pl^ is an easy observable to measure 
numerically, the numerical cost of this extra step is negligible. (However, since 0^^ is a global 
quantity, some extra work is needed when using parallel computer architectures.) 

From Eqs. (4^) and ( [4.7|) we see the main difficulty of the multicanonical method; we 



have to know the result we are after, Pcan., to some accuracy, before we can even start the 
multicanonical simulation! This requires some kind of bootstrap process, to be discussed 
below, in order to determine an initial guess for Pcan.; after the multicanonical simulation, 
we get an improved estimate for Pcan. from Eq. (|4.6| ). 

In the above discussion we implicitly assumed that the weight function W has been 
optimized for one particular value of m'jjj:. However, due to the factorization property of 
the m'jjj, term in the Hamiltonian, Eq. ( |2.5| ), we obtain the canonical probability Pcan. for 
a whole range of mj^rp values from a single multicanonical Monte Carlo run: for example, if 
the multicanonical weight function has been originally calculated with m'jprp = ml (and we 
have the resulting distribution Pmuca), we have 

Pcan.(m2; 4>l) oc Pn.uca(0 exp [ - W {c^D + ^{ml - ml)<pl] . (4.9) 

Naturally we can only determine Pcan.("^^; 0av) ^^^ values of (pl^ where our Pmuca determi- 
nation is accurate. Strictly speaking, when we perform a multicanonical run it does not 
correspond to any particular value of m^, since we can always absorb the m'^ term in the 
action into the weight function. 

How does the computational cost scale in a multicanonical simulation? Ideally, if we 
have guessed H^(0av) correctly, the system performs a random walk in the 0^^ range of 
interest, say, from 0^^^ to 0av2- Let us consider what happens if keep the range of (pl^ 
considered fixed as we increase the volume of the system. Now, if we require that the 
system "random walks" through the range a comparable number of times as the volume 
is increased, the computational cost is proportional to (0av2 ~ '/'avi)^^^- The factor V'^ 
appears because (pl^ is an intensive variable. For comparison, in a canonical simulation at a 
first order phase transition, the numerical cost rises as max02^(Pcan.(0av))/™^</>iv(-^can.(0av))- 
In a large box, as previously discussed, this scales as the exponential of an interface surface 
area, ~ exp(2cry^/^), which is a vastly more severe increase. 

In realistic situations there are often "hidden" barriers, which can hinder the random walk 
through the range of interest. For instance, in our case such a metastability occurs in large 
cubic volumes when the surface geometry changes: as shown in Fig. ^ there can be sphere ^-*- 
cylinder and cylinder ^^ slab transitions. At the transition point two different geometries 
have equal volume fractions and surface area. During a simulation the transition must occur 
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by the Monte Carlo algorithm finding a series of mixed phase geometries which smoothly 
interpolate between cylinder and sphere, which means that the surface area must increase for 
a fixed volume fraction. Thus, if we use extremely large volumes, the transitions between 
geometries become exponentially suppressed, even if the total multicanonical probability 
remains constant. (Of course for the study of bubble nucleation we will not have to study 
the range of 0^^ where such transitions occur; but when we determine Teq, or measure the 
surface tension (see below), it can become an issue. We also observe some metastability 
at the value of 0^^ where the dominant configuration changes from being homogeneous 
symmetric phase to a small broken phase bubble.) 

The success of the multicanonical method hinges on the accurate determination of the 
weight function. If we require that the resulting probability distribution Pmuca is constant 
up to factor of 2, say, then the weight function must be determined up to an accuracy of 
log(2) ^ 0.7. Worse accuracy in determining W significantly degrades the efficiency of the 
subsequent Monte Carlo, so such a requirement on the accuracy of W is actually necessary. 
The variation of W across the range of interest in this work is of order ?^ 100 (that is, we 
have to boost the probability of suppressed phase space regions by a factor of exp(lOO).) 
Thus, the weight function has to be determined to an overall accuracy better than 1%. 

We use a continuous, piecewise linear Ansatz for the weight function. We determine the 
weight function with an automatic iterative calculation procedure, using variations of the 
procedures presented in |]3^ , |30| . One approach is to choose a starting guess for W^(0av) (^^^ 



instance, a constant), and to perform a Monte Carlo under the (Markov chain) algorithm 
which would generate the distribution, Eq. ( [4.5|) . However, after each update sweep, W is 
decremented at the current value of 0^^. Thus, if some region of 0^^ is getting sampled 
very often, W is reduced there so it will be sampled less often. This procedure will cause 
W to evolve towards the correct form, but imperfectly because the recent history of the 
Monte Carlo is over-refiected in the resulting W. To fix this, the size of the decrements 
is reduced every time the evolution successfully explores the full range of (pl^ from bottom 
to top and back. When the total change to W, in the time the Monte Carlo evolution 
spans the range of (pl^, is negligible, the weight function has been determined with sufficient 
accuracy. Measurements of Pcan. then consist of two parts: first, we perform a run during 
which the weight function is iteratively improved to a required accuracy. Second, using 
this weight function, we perform a normal multicanonical run, which gives us the final 
probability distribution. The determination of W typically accounts for 30-50% of the total 
computational effort. 

C. An application: surface tension 

We illustrate how the multicanonical method works by measuring the surface tension 
(T - i.e. the free energy/area carried by the phase interface - with the histogram method 



45| . This has become the standard and well-understood method for computing the surface 
tension in a variety of lattice theories, including work closely related to ours, SU(2) gauge 
+ Higgs theories |]3|,|46|j9[] and effective theories for the MSSM ||34|| . 



As discussed in Sect. [Li A| , at the phase transition value of m^jj., where the symmetric 
and broken phases are equally probable, the mixed phase configurations with approximately 
equal volume fractions of symmetric and broken phases are exponentially suppressed (Fig. ^. 
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The suppression is proportional to exp(— Fsurfacc/^) = exp(— cr x Area/T). This is seen as a 
valley in the probability distribution of 0^^, see Fig. p. 

For the interface tension measurements it is advantageous to use lattices with cylindrical 
geometry, L^ ^ L^. = Ly. Because we use periodic boundary conditions, there will be at 
least two interfaces which span the lattice. The cylindrical geometry makes the interfaces 
tend to form parallel to the (x, ?/)-plane; and L^ should then be long enough so that the 
two interfaces do not interact appreciably. This is seen as a flat minimum in the probability 
distributions. 
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FIG. 6. Probability distributions of <i)%^ measured from cylindrical lattices at /? = 4:/{g^Ta) = 9 
and X/g'^ = 0.036. 



We perform a multicanonical Monte Carlo with the same lattice action and update as [Q, 
using the improved relation between lattice and continuum parameters found in 0. In Fig. |^ 
we show the probability distributions from 16^ x 128 and 24^ x 80 lattices at A/{g'^Ta) = 9 
and X/g"^ = 0.036,0 measured at the critical m^, which is determined by requiring that the 
symmetric and broken phases have equal probabilistic weight. (Hence, the technique also 
provides an accurate determination of the the critical value of m^, which we will need for 
comparison later on.) The larger lattice is not quite long enough to have the fiat central 
part the smaller lattice has. Note the striking difference in widths between the symmetric 
(left) and broken (right) phase peaks; this reflects the large ratio in the 0^^ variance in the 
two phases, discussed in subsection [IV A . 



^Corresponding to "unimproved" lattice Pg = 9.6674, x = X/g^ = 0.0389. For the lattice to 
continuum relations for m? and 0^ see fp|l. 
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FIG. 7. Geometric shape of a phase interface, from a Q2/g'^T across box. The height of the 
fluctuations goes as {T/a)'^''^; Fourier analyzing and averaging over hundreds of such surfaces can 
yield an accurate determination of a. 



Spacing a/g'^T 


Volumes used 


a 


4/7 


722 X gg^ gQ2 X ^20, 108^ X 144 


(0.0749 ± 0.0008)g^T^ 


4/9 


60^ X 144, 96^ X 144, 108^ x 160 


(0.0758 iO.OOM^T^ 


4/12 


120^ X 160, 1322 X 180 


(0.0733 ±0.0025)^^r3 



TABLE II. Surface tension as a function of lattice spacing. 

The surface tension is obtained from log[Pmax/-Pmm]/(2-^2) -^ a/T, as l^ ^ cxd. In 
practice, the infinite volume value of a is reached in such large volumes that finite volume 
analysis becomes necessary. Following Ref. |^ , ^ , we fit the data with the Ansatz which 



takes into account the translation modes of the surfaces and capillary fluctuations: 
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+ - log L^a — XogL^a + const. 



(4.10) 



The result of the fit from these two lattices is 



irp-i 



a = (0.079 ± 0.004)/T' 



(4.11) 



We also obtain the equilibrium m'^, and the difference in (pl^ between the two phases, which 
in physical units is Acf)^^ = ^^^(broken) — ^^^(symm) = 2.53g^T^. 

The surface tension can be obtained much more economically with an alternative method 
due to Moore and Turok |^|. This method is based on analyzing the spectrum of the 
transverse fluctuations of the phase interfaces; the magnitude of the fluctuations is inversely 
proportional to 



I a/T. We refer the reader to the reference for a complete discussion. We 
apply this method using multicanonical tools to sample conflgurations in a very large box, 
but now choosing W{(j)1^) to very strongly prefer 0^^ within 5% of the average between 
symmetric and broken values; hence the volume always contains large regions of each phase, 
with two approximately planar interfaces separating them. We show an example of such 
an interface in Fig. 0, and present the determined surface tension, as a function of lattice 
spacing, in Table |l[ The results agree within error with the histogram method. If we 
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extrapolate the values given in the table to a ^ assuming 0{a'^) errors, as should be the 
case since we use an 0(a) corrected lattice-continuum match, we obtain the result 



a = (0.0749 ± 0.0027) c/^T^ , a = , (4. 12) 

with a lattice spacing dependence which is small and consistent with zero. This indicates 
that our lattice spacing errors are under control. 



D. Results: probability distribution 

Let us now turn to our main problem, the determination of the probability distribution 
in the region relevant to bubble nucleation, see Fig. |^. The procedure is very similar to 
the surface tension calculation with the histogram method described above. However, there 
are two crucial points where it differs: 1) we need very large, preferably cubical volumes in 
order for the bubbles to fit in the lattice comfortably. The size makes it next to impossible 
to compute the full probability distribution from the symmetric to the broken phase with 
our computational resources. However, 2) we need the weight function only in the range 
0^y(symm) < (^^^ < (O.850^y(symm) + 0.150^^ (broken)), as discussed in Sect. |1IQ . This 



guarantees that we do not yet enter the "cylinder" and "slab" regions of the phase space, 
see Fig. §. Using the random walk argument, calculating the distribution in this restricted 
range only requires a factor of (0.15)^ ~ 0.02 of the resources needed for the full weight 
function.]^ 

In Fig. 1^ we show the probability distributions from a 124^ lattice, with 4/{g'^Ta) = 
9. The result of a single multicanonical run has been reweighted to 3 different values of 
mfjj^, namely (0.008, 0.009, 0.010)f7^T^ below the equilibrium value. The critical bubbles 
correspond to the minimum locations of the probability. Note that the value of (pl^ for the 
critical bubble moves to smaller values as we increase the supercooling; this is because the 
critical bubble gets smaller at larger supercooling, and 0^^ is a volume average. We also 
have results from a 92^ box at A/i^g^Ta) = 7, to study lattice spacing dependence, and from 
a 120^ box at 4/{g'^Ta) = 9, which gives weak information on volume dependence. The 
ratio of probabilities between the metastable minimum and the critical bubble, Pcan. (0av — 
0av.c)/-Pcan.(0av = 0av(syHmi)); IS plottcd for cach lattice as a function of 5m^ = m^^ — 
'>Tijjrp[equi\ih) in Fig. ^. This figure corresponds to the difference in height, in Fig. §, between 
the local maximum and local minimum. For each curve in the figure, the statistical error 
bars are about ±1, with strong correlation in the error along the curve. The finer spacing 
lattices agree within errors. This is a (weak) check on volume dependence, but it is also 
a check of the code, since the 120^ and 124^ volume computations were performed with 
completely independent sets of code, on machines of different architecture. The coarser 
lattice data differs by between 2% and 3%. This probably represents lattice spacing errors; 
if we extrapolate assuming a^ errors (the first order which should be present, due to our 
lattice improvement), we estimate the difference between the finer lattice and the continuum 



^'^This discounts the "barriers" at the bubble <-> cylinder <-^ slab transitions, which would make 
the full computation even more costly. 
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log(P) versus 0^^ at 3 supercoolings 
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FIG. 8. The probability distribution for 0^^ at three values of 6m'^, —O.OOSg^T^ (solid), 
-Omdg'^T'^ (dashed), and -O.OIO^^T^ (dotted), for 124^ lattice at g'^aT = 4/9. In each case 
the local maximum is the supercooled symmetric phase, the local minimum is the critical bubble. 
The three curves are obtained by reweighting the same multicanonical run. Greater supercooling 
leads to less suppression of the critical bubble. 
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FIG. 9. Log of ratio of probabilities logPcan.(nieta)/Pcan. (crit), as a function of 6m'^ the 
supercooling from the equilibrium m^. The solid line is the 124^ lattice, and the dashed line is the 
120^ lattice, at {A/g^Ta) = 9. Each has a statistical error of ±1, so they agree within expected 
error. The dotted line is the data from the 92^ lattice at (A/g'^Ta) = 7. Its disagreement represents 
definite, but small, lattice spacing error. 
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is 1.5 times as large as the difference in the curves. Where the ratio of probabihties is e^°°, 
this would be a correction of e'* in the nucleation rate, or about a 2.5% correction in 5m?. It 
remains to integrate the area in the symmetric phase minimum, to compute the dynamical 
contributions, and divide by the volume, to convert this result into the real time rate. 

E. Dynamical prefactor, tools and calculation 

To determine the real time rate for nucleations, we now have to perform real time evo- 
lution on each of a sample of configurations with 0^^ = ^^vc- ^^ S^^ ^^^ sample of such 
configurations by, first, choosing a 5m^ to consider, and, next, by performing a multicanon- 
ical Monte Carlo, as just described, and recording those configurations for which 0^^ lies 
within a narrow tolerance of (pl^ q (which depends on 5rn?, see Fig. g). In fact we can speed 
up the sampling process by choosing a weight function W{(t)\^) which favors (pl^ = ^^vC 
even more strongly than the one used to determine the probability distribution in the last 
section. Then, we must study the real time evolution of each configuration in the sample, 
at the thermal Higgs mass ml^ — 5m?. 

As discussed in subsec. pi B| , the appropriate dynamics, at leading log, are Langevin dy- 



namics. The gauge fields evolve according to Eq. ( |3.!J| ) , and the Higgs fields also evolve under 
Langevin dynamics, with a much (parametrically) faster time scale. In continuum notation, 
this means evolving the fields under the following Langevin field equations (normalizing E 
so E"^ j2g^ appears in the action), 

^ei^r(^) = -^]|(^^+er(^,t), (4.13) 

d 

aeiA$(a:) = -^^^t(^^ + ^*(^' ^) > 

(C(x, t)e)(x', t')) = 2aeir5,,(5"^(5(x - x')5{t - t') , 
(e*(a;, tKiix', t')) = 2r]a,,T 15{x - x')5{t - t') , 



where H is given in Eq. ( p.7|) , 1 is the identity in component space for the Higgs field and 
T] is the ratio of the speeds of Langevin evolution, which is parametrically ~ l/a^, and so 
should be taken large. 

It is possible to perform numerical Langevin evolution on lattice fields, but it is slow 
and unnecessary; any dissipative update will do, if the relation between the number of 
updates and Langevin time is known. Hence we use the heat bath algorithm to update the 
gauge fields; the relation between the number of heat bath updates, and the time scale t 
in the equations above, is discussed at some length in |TB[. At leading order in small a the 



relation is that, for random order heat bath updates of the lattice links, n updates per link 
corresponds to At = a^o"ein/4. (Note that this relation is specifically for our choice of lattice 
action, namely "Wilson glue" ; it would be different for an improved action. It can also differs 
if the sites are updated in a specific, rather than random, order, and in fact depends on the 
order of update.) We update the Higgs fields with a mixture of the over- relaxation algorithm 
presented in Q] and a Higgs field heat bath algorithm. Note that our "real time" evolution 
algorithm can be viewed as a canonical Monte Carlo evolution algorithm at m? = ml —5m'^] 
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FIG. 10. A tunnelling trajectory and a trajectory which does not lead to tunnelling, measured 
from 120^ volume and at supercooling 6m? = — 0.0082 g^T^. The horizontal dashed line is the 
critical bubble value of (/)^v Two half-trajectories (positive and negative timestep values) are 
evaluated starting from a configuration at the critical value of 0^v5 ^^^ glued together at timestep 
to form a full trajectory. 

hence there is no concern that it somehow spoils the thermodynamics (as might happen for 
Langevin or Hamiltonian evolution with a finite time step, due to time step size errors). 
As described in subsection |11B|, we calculate the "dynamical prefactor" d, Eq. (|2.8|), by 



evolving a critical bubble configuration both forward and backwards in time, long enough 
to see whether the system evolves either towards the symmetric or the broken phase. Since 
forward and backwards Langevin evolution are equivalent, in practice we generate a few 
Langevin trajectories from each initial critical bubble configuration; each pair can be joined 
together to form a full trajectory, see Fig. [l^, and by considering all pairings we somewhat 
improve the statistics. The dynamical prefactor is then the expectation value 



N, 
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6, 



tunnel 



t^'^j- traj 



^ # of crossings 



(4.14) 



where (5tunnei is 1 if the trajectory leads to tunnelling, otherwise, and (# of crossings) is the 
number of times the trajectory crosses the critical bubble value of 0^^. When computing the 
error in the determination of d we must account for the dependence of the several trajectories 
involving a common critical bubble. 

One legitimate concern is that the Langevin dynamics we consider are only correct at 
leading log, and as argued in 
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the treatment should break down when there is a large 
Higgs field condensate. Hence our treatment of the dynamics may not be better than an 
0(1) treatment. However, our determination of the thermodynamic likelihood of a critical 
bubble was only good to ±1 in the exponent, so an 0(1) error in the dynamics is no worse; 
in any case, since the nucleation rate is ~ e~^°°T^, a factor of 2 error in its determination 
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only represents a 1% error in the exponent, and somewhat less than a 1% error in the 
determination of 5rn?, as seen in Fig. p. It would be possible to do a better job by using the 
full HTL dynamics by the technique developed in [jl^; however, this approach is much more 



numerically expensive. It also requires the inclusion of the Aq field in the thermodynamic 
treatment, and it could be difficult to eliminate finite time step errors in the dynamics, 
which could make the thermodynamics explored by the evolution slightly incorrect. 

The reader might also be concerned that the dynamics of a broken phase bubble will not 
have a good limit as rj is taken large. If the Higgs fields are allowed to evolve much faster than 
the gauge fields, won't the bubble either collapse or expand on a time scale set by the rate 
of the Higgs field evolution? The answer is, no. If we choose a starting configuration with 
a broken phase bubble, and we evolve the Higgs fields without allowing the gauge fields to 
evolve, the bubble does not collapse but remains a critical bubble indefinitely. It is essential 
that the electroweak phase transition is a radiatively induced phase transition; the state of 
the gauge field fiuctuations, alone, are sufficient to indicate which phase a configuration is 
in, and the Higgs condensate cannot expand into the symmetric phase, or collapse inside the 
bubble, without the gauge field fiuctuations changing as well. To demonstrate this point, we 
have performed an 77 = cxd evolution, meaning an evolution in which only the Higgs fields, 
and not the gauge fields, are allowed to evolve. It is compared to an evolution in which both 



evolve, but the gauge fields evolve much more slowly, in Fig. 11. Naturally, the value of 



av 

changes during each evolution; but in the evolution with frozen gauge fields, it just bounces 
around a central value, and the bubble is stable. We should expect, then, that (|A0^y/At|) 
and d will depend on 77. But, as we must check, the product should have a finite limit. We 
would also like to know how the product (|A0^y/At|) d depends on drn?] we expect that the 
dependence is weak. 

To check this we have used two values oi t], t] = 5 and rj = 10, and measured d and 
{\A(j)l^/At\) for each (on a 124^ lattice at 6171^ = —0.0082g'^T'^). The results are presented in 
Table |T^, which shows that the product d(| A0^y/At|) is independent of rj, within numerical 
errors. We have also re-analyzed the same set of trajectories, sampling (p^^ half as often; we 
find as expected that d and (|A</)^y/At|) each depend strongly on the sampling rate, but 
the product does not. 

We have also varied the degree of supercooling, to confirm that the dependence of the 
dynamical prefactors is not very strong. In the 120^ box we have results at 3 different 
values of supercooling: 5rn? / {g'^T'^) = —0.0082, —0.0093 and —0.0104. These correspond to 
significantly different bubble probabilities, as shown in Table |V[ the largest supercooling 
corresponds to bubbles which are ~ e^^ times more likely to nucleate than with the smallest 
supercooling. The bubble probability density, Eq. (|2.7]), is readily evaluated by integrating 
the data shown in Fig. §. However, as expected, the dynamical factors are seen to be fairly 
stable throughout this region, varying only by approximately a factor of two or three, which 
is in practice insignificant when compared with the differences in probability. 

Note the units on (|A0^^/At|). The units on the first term in Eq. (|2.6| ) are the same 
as l/0av ~ ^/c^wT^] while to get a rate per unit volume we must divide by the volume, 
1/V ~ chIjT^. Hence the parametric appearance of the nucleation rate is oc a^T^log(l/5f), 
using (Tel ~ T/\og{l/g). This arises simply from the relation between dimensionless lattice 
quantities and physical quantities. 
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FIG. 11. Two time histories in the 124^ volume, both starting at the critical bubble: that 
in which the gauge field is not updated (blue, stays nearly constant) and that in which the gauge 
field is updated (black). When both fields are updated, the critical bubble can grow, or, in this 
case, collapse; if only the Higgs fields are updated it does not evolve. 
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TABLE III. Dynamical information (|A0^y/At|) and d, varying rj and sampling each data 
series at two rates (all on a 124^ lattice at Sm'^ = — 0.0082g^T^). Each pair of data at fixed 
r/ comes from the same set of trajectories; comparison shows that d and (|A0^y/At|) depend 
strongly on sampling frequency, but the product does not. The (fully independent) data sets with 
different rj agree within error for (|A0^y/Af |) d, showing this quantity has good large rj behavior. 
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TABLE IV. The nucleation rate calculated with 3 different supercoolings, in a 120^, g^Ta = 4/9 
lattice. Pc is the probability density that (f)^^ = (t>%j q, calculated from Eq. (2.7). The nucleation 
rate in the fifth column is Pc(|AC/^*l>d/2. 
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FIG. 12. S, as defined in the text, as a function of 6m'^, the degree of supercoohng, for 
the g^Ta = 4/7 (dotted) and g'^Ta = 4/9, 124^ lattice (sohd) data, each using {\A(j)l^/At\)d 
as evaluated at 6m'^ = —0.0082g^T^. The points are the 120^ data, using the three evaluations 
of (|A(/)^y/At|)d at the values of dm? where they were evaluated. Error bars, not shown, are 
dominated by the error in the determined probability distribution, and are about =bl in S. 

F. Numerical results and relation to cosmology 



The previous two subsections contain all the ingredients needed to determine the real 
time rate for bubble nucleation. It remains, first, to put the ingredients together, and second, 
to determine what value for the nucleation rate is interesting cosmologically. We write the 
bubble nucleation rate as 



rate 



2t^2' 



g'T 



m 



\ogil/g)alT'expi-S). 



(4.15) 



D 



For the realistic standard model values mj^ = {ll/6)g'^T'^ and a^, — 1/30, minus the log of 
the term in front (evaluating \og{l/g) using Eq. p.9|)) is 16, so the rate is exp(— (5* + 16))T^. 
Our result for 5* is shown, as a function of 6m^, in Figure |1^. This represents our final 
numerical result. It is unfortunate that the numerical effort is too large to perform the 
calculation for several couplings, and we have also not considered a realistic set of parameters 
in the MSSM. 



What are the errors of 5* in Figure [12| ? Using the numbers in Table |^, for example, we 
see that the final statistical errors are strongly dominated by the errors of the probability 
distribution P. In addition, there is the systematical error of the real-time update evolution 
IV E| ). Both of these error sources are easily included by a ±2 error band around 



see sect. 



S in Figure [T^. 

We want to know at what value of S the phase transition occurs, so we can determine 
Sm"^ and therefore the amount of supercooling. The relevant picture is discussed in |Q . The 
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bubbles of broken phase which convert most of the volume into the broken phase, nucleate 
over a characteristic period of time tnuc, and with a mean separation d^uc, which is also the 
diameter they grow to. For one such bubble to nucleate per d^^^ volume in time tnuc, the 
nucleation rate must be ~ l/(tnuc'^nuc)- The time scale is set by how fast the nucleation 
rate is changing; tnuc — (dS/dt)^^ ~ {H(T)TdS/dT)^^, with H{T) the Bubble's constant 
at temperature T. The separation is d^uc ~ ^s^nuc, with Vs the sound speed ~ 1; this is 
because a bubble is preceded by a shock front propagating at approximately Vg, which heats 
the plasma, and nucleations are suppressed in the heated region. Hence the nucleations take 
place when 

rate^(^i7^) . (4.16) 

From Fig. ^, and from the mass-temperature relation, Eq. ( ^.31 ), we see TdS/dT ~ 2 x 
10^; while from the Friedman equation, H = J4:Tr'^g^/4:5{T'^/mpi) ~ e~^^'^T. Hence the 
interesting value of S" is 16 + S* ~ 106, or 5 ~ 90. The g'^Ta = 4/7 data show this 
value is obtained at Sm^ = (—0.00880 ± 0.00010)g^T^, while the finer lattice data give 
(— 0.00864 ±0.00010)5'^T^. An extrapolation to zero lattice spacing, assuming Ola^) errors, 
gives Sm'^ = —0.00840 ± 0.00020. Again, this is for X/g'^ = 0.036; naturally the amount of 
supercooling is strongly dependent on the strength of the phase transition, and hence on 
X/g"^. This is our final numerical result. 

V. OTHER APPROACHES 

Since we have only computed the nucleation rate nonperturbatively for a single set of 
parameters, the main application of this work is as a benchmark for studying the performance 
of other, less first principles means of determining the nucleation temperature, which have 
previously been used in the literature. For this purpose, we will apply a few of these 
techniques to the current set of parameters, to see how accurately they determine Sm"^. 
(This actually gives an optimistic appraisal; the dependence of S on dm"^ is strong, so the 
relative error in S the log of the nucleation rate is typically almost twice the relative error 
in Sm'^.) 

A. Thin wall approximation estimate 

The idea of the thin wall approximation is to hope or assume that the critical bubble is 
"thin walled," which technically means we assume 3 things: 

1. the bubble radius is much greater than the thickness of the phase interface, so we can 
treat the interface as an infinitely thin geometric surface; 

2. the surface tension of the interface, at the nucleation temperature, is the same as at 
the equilibrium temperature; 
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3. the free energy difference between the phases is given by the leading order expression 
AV^ = lAT/T, with / the latent heat at equilibrium. This neglects the change in I 
between the nucleation and equilibrium temperatures. 

In the limit of small supercooling, all three approximations become exact. Hence the thin 
wall approximation is appropriate if we are interested in very small 6m'^, where the critical 
bubble free energy is huge. The approximation is in expecting it to continue to work down 
to where the critical bubble free energy is S* ~ 90T. 

The thin wall approximation is that a bubble of radius r will have energy 

E{r) = 47rrV - —r^AV , (5.1) 

o 

with (7 the equilibrium surface tension and 

We find the extremum over all r, 

dE 2a levTCT^ 

We now estimate that the nucleation rate will be 

rate = exp(-E/T) f ^^ al \og{l/ g)T^ , (5.4) 

where the technique really provides us no way of getting the non-exponential term we write, 
but we guess that this is correct on parametric grounds. 

For comparison with the case we have studied numerically, we solve for the value of bvn? 
required to make EjT = S* = 90, using the nonperturbatively determined values of a and 
A02^(Tcq) from the last section, o = 0.075g^T^ and A02^(Teq) = 2.53g^T^. Substituting 
into 

gives 5m^ = — 0.007 Og'^T'^. This is low by about 20% from the nonperturbative value. 

Alternately, we could ask how accurately the thin wall approximation predicts the nucle- 
ation rate at the actual value of drn^. Plugging in the value of (5m^ found in the last section 
gives S = 62, which is off by almost a third. 

It is clear that the thin wall approximation is doomed to errors of this magnitude when 
applied where E/T ~ 90. Note that the radius of the critical bubble, for 6771"^ determined 
above, is 



&'■ 



But it is impossible for the thickness of the interface to be less than ~ JT/a; even if 
the intrinsic thickness were somehow thinner, the surface fluctuations would generate such 
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a thickness. Hence r never greatly exceeds the surface thickness if we are interested in 
E /T ~ 90. However, all considering, the thin wall approximation (given nonperturbative 
inputs) does pretty well. 

We should comment that, although the thin wall approximation requires nonperturbative 
inputs, it is much easier to determine A0^^ and a on the lattice than to directly compute 
the nucleation rate. In particular, we could think realistically of doing a "scan" of A^^^ 
and a at several parameter values, in the standard model or one of its extensions, but to 
do the nucleation rate at numerous values of parameters would probably be beyond what is 
currently a reasonable amount of numerical effort. Hence the thin wall approximation may 
be a reasonable approach to getting "rough and ready" nucleation information. 

B. Perturbative estimate 

The traditional method for determining the bubble nucleation rate in the context of a 
perturbative treatment of the strength of the phase transition (see for instance pO|-p^) is 
to approximate the free energy (effective action) to be the tree kinetic term for the (gauge 
fixed) Higgs condensate 0, plus an effective potential term computed at some order in the 
loop expansion, 

E = jd'x{^-{Vct>f + V{ct>)) . (5.7) 

It is also possible to consider radiative corrections to the Higgs field kinetic energy, as we 
discuss in the next subsection.|^ 

We begin with the simplest possible estimate. We take the one loop effective potential, 
Eq. ( p.8|) , reproduced here for convenience: 

^iioop = ^0^-^0^T+V. (5.8) 

/ ibvr 4 

Here we have neglected scalar loops, which give a contribution down relative to the (f)^ term 
by (A/^f^)^/^. We can see from Eq. ( pISD that at 



^^We note that a full-fledged analytical computation of the nucleation rate in the spirit of the 
Langer method is very difficult in radiatively induced first order transitions. The problem is how to 
distinguish the fluctuations which give you the effective potential from the fluctuations of the bubble 
in this potential. The results from a cubic anisotropy model (a simple spin model) calculation by 
Strumia and Tetradis |5^ display a dramatic dependence on how the separation of the fluctuations 
is done; since the physics cannot depend on the separation of fluctuations, this effect must be an 
artifact of the approximations done in the calculation. On the other hand, if there is a sufficiently 
strong first order transition already at tree level, the Langer method is relatively straightforward 
to apply reliably |54-56|. 
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there are two degenerate minima with 

(f) = and (f) = (j)Jl loop) = ^ . (5.10) 

SttA 

At best perturbation theory is an expansion in ~ (g'^T/mw) ~ qT / ip, hence in \/ g^. Two 
loop corrections will give a correction to 0o of order {X/g'^)(j)Q, which is larger than the scalar 
loop contribution we neglected above. This justifies neglect of the scalar loop.Q 

One then assumes, reasonably, that the critical bubble, the saddle point of Eq. ( |5.7| ), will 
have spherical symmetry, and looks for the saddle point of 

oo /I 

"2 1^ At^W'i 



E = ^7tJ^ r'^-{dr(j){r)y + V{r {<!))) j (5.11) 

over all 0(r), with the boundary condition that 0(r) go to zero at large r, so we have a 
bubble in the symmetric phase. One finds the saddle point free energy as a function of 
rnjjrp and looks for where it takes on the desired value, say E/T = S = 90. The difference 
from the equilibrium m'j^rp is the degree of supercooling we are after. (As in the previous 
section we will simply assume that the bubble zero modes etc. give a dimensionful prefactor 
of ~ a^T^. Since even a change in the numerical value of this prefactor by a factor of 1000 
would represent only a change of 7 in the log of the nucleation rate, our results depend only 
weakly on this treatment.) 

We are not aware of a purely analytic way to find the saddle point of Eq. ( |5.11| ), but 
the well known overshoot-undershoot algorithm is both efficient and accurate. We find, for 
X/g^ = 0.036, that 



5m' = -0.0058(?^T2 , aeq^mbnum = r V2Vd(f) = ,,,/J,,,,, = 0.0302^^T3 




3072v^7r3A5/2 



(5.12) 



In comparison to the nonperturbative values these are both fairly far off. The degree of 
supercooling is 2/3 of the right value and the surface tension is 40% of the nonperturbative 
value. One loop perturbation theory is NOT accurate even at X/g' = 0.036, though it is 
not completely wrong. 



^^As an aside we mention that the expressions usually written for scalar loops (refs. |57| , |58| , 
though not ref. ||5^) must be inconsistent because they depend on mj^rp in a way which does not 



satisfy Eq. |2.5| . The problem is that they were derived using a Higgs mass in loops which does not 
correspond to the curvature of the potential; the value of rriff used in loops includes the second 
derivative of the Xcp^ potential term but not of the negative cubic term, which is of the same size 
in the interesting region and cannot be neglected in any consistent approximation scheme. If we 
include terms in the effective potential in the spirit of an expansion in X/g' then no reference to the 
scalar self-coupling or m'jjrp appears in any loop induced term until 0((A/gf^)'^'^), when an infinite 
class of diagrams must be resummed, corresponding to a single Higgs loop but including iterated 
one loop mass corrections from gauge boson loops within the 3D theory (not just the nonzero 
Matsubara frequency loops). 
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Potential used 


Wave Function Z 


-5mV(7^T2 


ajg'T^ 


1 Loop 


Tree 


0.0058 


0.0302 


2 Loop 


Tree 

■^exp 
■^pade 


0.0168 
0.0109 
0.0115 


0.088 
0.072 
0.074 


"nonperturbative" 


Tree 

■^exp 
■^pado 


0.0151 
0.0110 
0.0114 


0.097 
0.083 
0.084 


Nonpert. Result 




0.0084 


0.075 



TABLE V. Supercooling and surface tension, using several analytical or semi-analytical ap- 
proaches, compared with the full nonperturbative answer (last line). The meanings of Zpade and 
Zexp are explained in the next subsection. The one loop result gives too small results; all other 
results give too large an answer. 

We therefore go on to the two loop effective potential. Neglecting powers of A and setting 



vTiHT -c mvF to zero within loops, the new term in the effective potential is |60 



V, 



2 loop 



Vi 



51 



loop 



5127r 



.^V'T^ log(</./^T) 



(5.13) 



where the choice of gT inside the logarithm is a convenient choice of renormalization point; 
a different choice can be absorbed into a shift in m'jjrp. The justification for neglecting scalar 
effects is the same as before; for more discussion see footnote |^. 

Including this term in the effective potential, we find that, while 0g = 2.19 is still 
smaller than the nonperturbative A0^^, we now get too large a surface tension and too 
much supercooling, by a substantial margin: 



4:rr3 



cr2 loops = O.OSSg'^T 



Arr2 



6m\2 loops) = -0.0168^^T 



(5.14) 



The amount of supercooling is too large by a factor of 2. At the value of 5m? found 
nonperturbatively, 5*2 loop > 260 is almost 3 times too large. All the results for various 
effective potentials are summarized in Table 0. 

It is also possible to define a nonperturbative "effective potential for <^," as follows. We 
first find the largest volume for which the configurations with 0^^ intermediate between 
the pure phase values, do not show phase segregation. This turns out to be a surprisingly 
large volume. To see this, recall why, in a larger box, a configuration with 0^^ intermediate 
between homogeneous phase values is a mixed phase configuration, rather than an extensive 
region with the intermediate field value. The effective potential we will eventually derive is 



shown in Fig. |T^. The state half way between minima has negative curvature. Infrared Higgs 
field fluctuations are therefore spinodally unstable. If we force 0^^ to maintain its value, 
this prevents the zero mode from growing, ie. it keeps from shifting value homogeneously 
through the box; but small k modes will be unstable. Any mode with k < uj_, with 
uj_ = \/ —V" evaluated at the unstable point in question, will be unstable to grow. However, 
in a finite volume, there is a discrete spectrum of k modes. The lowest nontrivial mode has 
k = 2tt/L] so in any volume smaller than L = 27r/u)-, configurations with intermediate (pl^ 
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0/ gT 



FIG. 13. Left: free energy (— In (probability)) distribution as a function of (/)^-y in a 40^ box 
at a = 4:/9g'^T (/3 = 9). The numerical errors are < 0.2 and not shown on the plot. Right: plot 
at left interpreted as an effective potential for cp, as described in the text. The dotted line is the 2 
loop perturbative result, included for comparison. 



will be homogeneous. In larger boxes they will be inhomogeneous, containing a region closer 
to one phase and a region closer to the other. 

Now, we determine F(0^y), meaning — log(Prob(0^y)), in such a "large but not very 
large" volume. In practice we use a 40^ lattice at /3 = 9, meaning a volume 17.8/g'^T on a 
side. This is pretty big, though the volumes we used to study critical bubbles were typically 

F/V and = 



3 times longer on a side. We then write V 



9 



^U 



s phase), with 



U 



s phase) the lower local minimum of F{(j)l^) at Teq. That is, we throw away the part 
of the picture on the left in Fig. |l^ which lies at smaller 0^^ than the first minimum, and 
rescale the x axis for the rest. We take the result to be a "nonperturbative measured V{(f))" , 
shown at the right in the figure. 

The advantages of getting an effective potential in this way are that it should automat- 
ically get the right latent heat, and will show the disappearance of the phase transition at 
large X/g"^. The disadvantages are that it is numerically expensive (though less so than 
a direct determination of the nucleation rate), is somewhat arbitrary (how exactly do we 
choose a volume, for instance?), and is still only approximative. In particular, it is not at all 
clear that it is reasonable to assume that as defined above will have canonically normalized 
gradient energies. 

The "nonperturbative effective potential" has a larger separation between minima, but 
the height of the barrier rises relative to the two loop potential more weakly than as 0o- ^^ 
a result it gives a slightly lower degree of supercooling, Sm"^ = —0.0151, as summarized in 
Table 0. Note however that the surface tension is further off than in two loop perturbation 
theory. What this exercise teaches us is that the difference between the perturbative and 
the nonperturbative values of a and 5m^ are not primarily because of the limitations of the 
2 loop effective potential. 



38 





FIG. 14. Diagrams which lead to 0{X/g^) important Higgs wave function corrections. The 
vertices in the diagram at right only exist because there is a condensate. 

C. Perturbative estimates including wave function corrections 

We have just found that estimates of the degree of supercooling and of the surface tension 
were not substantially improved by passing from the one loop to the two loop effective 
potential, using the tree level Higgs field kinetic term in each case. In fact this result is 
not surprising. While the two loop effective potential can tell us the latent heat / oc 0^^ at 
next to leading order in \/ g"^., the procedure we used is still correct only at leading order 
in \/ g^ for determining the quantities we want, because it does not incorporate thermal 
corrections to the Higgs field wave function. In this subsection we see what happens when 
we incorporate Higgs wave function corrections which make the calculation complete to next 
to leading order in \/ g'^. 

Using the one loop effective potential, we can get the parametric estimate 0o ~ g^T/X, 
so mw{4') ~ g^T/X in the broken phase and inside the wall. The wall thickness, on the 
other hand, is set by the Higgs boson mass in either phase, L„aii ~ ^/rriHT ~ g^T/\/X. 
Therefore, up to a correction suppressed by a half power of X/g^, we may take the W bosons 
to be heavy compared with the reciprocal wall width, and treat the Higgs condensate as a 
homogeneous background for them. For this reason it is possible to incorporate the leading 
non-potential correction as a wave function correction for the Higgs field, 

IW? => ^W?. (5.15) 

Expanding the diagrams shown in Fig. O to second order in external momenta (treated as 
much less than mw), we find a correction to the Higgs gradient term which is, in Landau 
gauge, 

Z{<P) = 1 + 6Z{<j>) = 1 - ^^ = 1 - il^ , (5.16) 



in agreement with the expression found by Bodeker et. al. JUT]] (see also [^, where these 
wave function corrections have been considered at length). 

Both the corrections from Higgs fields, and from higher order in the p -C m\Y expansion, 
will give at most 0{X^^'^/g^) corrections to a and 5m'j^rp; whereas we see from the estimate 
for 00 that the above term contributes at 0{X/ g'^). It will drive down both a and (5m^, since 
it lowers the gradient energy contribution to the energy, and so makes it cheaper to have 
bubbles or interfaces. We should mention that to compute the {X/ g'^Y^'^ corrections, we 
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would have to perform a fluctuation determinant, since the Higgs field mass and the width 
of the interface are parametrically the same. 13 

Unfortunately Z{(j)) goes to zero at a finite value of 0. At the value of where this 
happens, the condensate is so weak that perturbation theory is also breaking down; the 
failure signals that neglected higher order effects become essential. To deal with this, we 
will make an Ansatzioi those effects, chosen to maintain the correct large behavior of 2'(0) 
but to prevent Z{(f)) from going to zero. We have considered two choices; to approximate 

Zexp = exp{-5Z) , and Zp^dc = f^ • (5.17) 

i — oZ 

Since we do not know the higher order behavior of Z, we do not know which of these is more 
sensible. If the answers we get depend strongly on which one we take, that is an indication 
that the perturbative expansion is failing and we cannot trust either. 

When we re-compute dni^ and a, including these wave function corrections, we get a 
result which is much closer to the nonperturbative value, see Table |V|. Further, the choice of 
how to resum higher terms in Z{(j)) appears not to matter much. However, the supercooling is 
still over-estimated by about 25%. At 6171'^ where the nonperturbative calculation shows the 
phase transition takes place, where ^nonpcrt = 90, the 2 loop potential with Fade resummed 
self-energy corrections gives S = 143, more than 50% high. 

In summary, using a perturbatively computed effective potential but tree Higgs kinetic 
terms appears to work badly. Higgs wave function corrections are important and improve 
the performance of the perturbative treatment; they should be included in any subsequent 
work which tries to study electroweak bubble nucleation by perturbative means. However, 
even with wave function corrections, perturbation theory is still not a very accurate way to 
treat bubble nucleation. 



VI. CONCLUSIONS 

In this paper we have presented a technique for determining bubble nucleation rates 
in theories with classical infrared thermodynamics and dynamics, in a fully nonperturba- 
tive way-even where the rate is exponentially small. The technique can be considered a 
generalization of hanger's formalism |2^, which replaces the approximate saddle point ex- 
pansion of that method with a nonperturbative Monte Carlo calculation, and takes care to 
treat correctly the microscopic dynamics of the nucleation process. The procedure uses an 
interesting mixture of multicanonical Monte Carlo and real time techniques. Within the 
context of the dimensional reduction approximation for the thermodynamics, and Bodeker's 



^^Such a fluctuation determinant calculation has been performed by Baacke [^, and has also 
recently been considered by Parnachev and Yaffe ^^ ; but since neither reference uses the two loop 
effective potential, their results cannot be considered correct through to 0{\/g^). It is not clear 
to us how to simultaneously perform the fluctuation determinant and to incorporate the two loop 
gauge contributions to the effective potential. Without resolving this question, the calculation of 
the fluctuation determinant is not justified in the sense of an expansion in X/g^. 
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effective theory for the dynamics |T^, our treatment is exact up to small and controllable 



numerical errors. It is useful when the microscopic physics is known well enough, and suit- 
ably amenable to a lattice treatment, to permit an accurate first principles Monte Carlo 
calculation. In particular, it can be applied at the electroweak phase transition. 

We have applied the technique to the electroweak phase transition in the standard model, 
at a value of the coupling which is just enough to preserve any baryon number after the tran- 
sition; the ratio of (3-D effective theory couplings) was taken as \/ g"^ = 0.036. The degree 
of supercooling is such that the thermal Higgs mass squared falls by 6m'^ = — 0.00840f7^T^ 
from its equilibrium value. 

The main value of this measurement is that we can compare it to the results of more 
traditional and less first principles calculations. We find that the most common method 
in the recent literature, using the two loop effective potential but tree level Higgs kinetic 
term, is very unreliable. For the parameters considered it overestimates the amount of 
supercoohng by a factor of 2, even though it gets 00; the length of the Higgs field in the 
broken phase, to within 10% error. If one considers the action of the critical bubble, it 
is even further off. Including Higgs field wave function corrections substantially improves 
the accuracy; the amount of supercooling is then only overestimated by about 25%. The 
remaining discrepancy cannot be attributed to the inaccuracy of the effective potential; a 
nonperturbative "effective potential for 0" is off by the same amount. On the other hand, 
a thin wall estimate, with nonperturbative inputs, does quite well-it is off by 20%, in the 
opposite direction. 

It should be straightforward to apply our technique to the more phenomenologically 
interesting MSSM or NMSSM, which can support viable baryogenesis. 
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APPENDIX A: DISCUSSION: INDEPENDENCE ON CHOICE OF 

MEASURABLE 

In this appendix we justify why the procedure presented in Section ^ of the main text 
"works," and in particular why the determined rate should be independent of the choice of 
measurable. 

Fix attention to a particular volume, and a particular value of T < Tgq, at which the 
nucleation rate is small. We must begin by clarifying the definition of the nucleation rate. 
First, we must be able to say whether or not a configuration is in the metastable phase. 
This requires that we possess some measurable in terms of which the free energy shows a 
"two well" structure, like Fig. |, with the probability at (C) exponentially smaller than at 
(A). For simplicity of notation let us assume that 0^^ serves this purpose. 

We can now draw two lines, roughly at {B) and (D) in Fig. |^; all we require is that 
they be well between (C) and the local minima, such that the likelihood to be at (B) is 
exponentially greater than the likelihood to be at (C), but exponentially less likely than to 
be at (A) (and similarly {D) is exponentially more likely than (C) but less likely than the 
broken phase). We say that a configuration is "definitely in the metastable phase" if it has 
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a value of 0^^ to the left of (B) , and is "definitely in the broken phase" if it has a value to 
the right of (D). 

The intuitive meaning of the nucleation rate is the following. Take the canonical ensem- 
ble, and throw out all the configurations to the right of (B), leaving only the ones which 
are clearly in the metastable minimum. Now carry out the time evolution for a "medium" 
amount of time, tmedium, exponentially longer than any microscopic time scale, but exponen- 
tially shorter than the time it takes for most of the metastable configurations to escape to 
the broken phase. At the end of the time evolution, look to see how much of the ensemble 
is definitely in the broken phase, ie. to the right of {D). Define that fraction, divided by 
^^medium, to bc the spacetimc rate of nucleations. 

This intuitive meaning of the nucleation rate will be our definition. Note that it is 
equivalent to the following. Consider the set of all trajectories of length tmedium, starting 
from any configuration (symmetric, broken, or in between) with appropriate Boltzmann 
weight. The nucleation rate per unit volume is 

P(symm —>■ broken) 

rate = — — — ^ r- , l-^-'-j 

»/imediumP(synim -^ symm) 

where (P) means the fraction of the trajectories satisfying the given condition. In other 
words, find what fraction of trajectories start in the symmetric phase and end in the broken 
one; and divide by the number which start and end in the symmetric phase, the volume, 
and the time. (The denominator should read P(symm — > left of (-D)), but since configu- 
rations between (B) and (D) are exponentially rare the difference is exponentially small.) 
This definition of a rate depends on our exact choices for (B), (D), and tmedium, but the 
sensitivity should be exponentially weak — unless there is some additional long time scale 
in the problem, besides the nucleation rate, in which case our technique probably fails. We 
can think of the nucleation problem in terms of sampling over the space of real time trajec- 
tories, and the sampling may be taken using the equilibrium ensemble. The goal is now to 
show that our technique correctly determines the number of crossing trajectories, relative 
to all trajectories which remain near the symmetric phase; and that this does not depend 
sensitively on our choice of measurable. 

From here on we will consider our technique applied to two classes of dynamics. The 
first is Hamiltonian dynamics, for a system where the phase space is the tangent bundle 
over the space of configurations. We further assume that the Hamiltonian is a function of 
space, plus (a constant times) the inner product function on the tangent space. In other 
words, the momenta appear quadratically in the Hamiltonian, and in any orthonormal basis 
for the momenta the quadratic form is a multiple of the identity. Probably these conditions 
are too strong, and we could work with any Hamiltonian system which could be written as 
a fiber bundle over a configuration space; but the treatment would be more complicated. 
The other class of dynamics we consider is Langevin dynamics, where the noise term is 
Gaussian, white, and additive. White means there are no unequal time correlators in the 
noise, in which case it is described by a weight function on the tangent space; we take that 
weight function to be of the same form as the momentum distribution just described for the 
Hamiltonian case. Probably one could extend our technique to the case where the noise is 
stronger in some coordinate directions, or has an amplitude which varies in space. However, 
such Langevin systems are notoriously subtle, see eg. |l65| , |66| . 
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Space of configurations 



ines of constant order parameter 
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Surface D 



This line of constant 

order parameter is the 

separatrix for that order parameter 

FIG. 15. A cartoon of how configuration space is foliated by a measurable-here the hnes of 
constant measurable do not coincide with the ones which determine the surfaces B and D. 

Consider two measurables, call them Oi and O2. (They could for instance be (pl^ and 
some other volume average of a local observable; or O2 could for instance be the maximum 
over all centerpoints for spheres of radius r, of /gp[^gj.g 2<l>''^$.) For each observable we can 
make a free energy plot and identify some least likely value of O, Ocrit- The measurable O 
is a map from the space of configurations to the real numbers; the subspace which gives a 
particular value of (9 is a codimension 1 surface in the space of configurations. Hence defining 
a measurable also defines a foliation of the configuration space. We will call the special 
surface, consisting of all configurations with O = Oait, the separatrix ior the measurable O. 
This notation is in keeping with [^], where many of the root ideas of our algorithm can be 
found. One would like the separatrix to separate the configurations more likely, under time 
evolution, to go to the symmetric phase, from those more likely to go to the broken phase. 

Note that [B) and {D) also give codimension 1 surfaces in the configuration space, 
namely the surfaces where 0^^ takes on the values (B) and (D); we will call these surfaces 
B and D. To one side of the surface B is the metastable phase, to the far side of surface 
D is the broken phase, and in between are intermediate configurations. For the measurable 
O to be useful, we require that the O separatrix carry all but an exponentially small part 
of its weight between the surfaces B and D. We give a cartoon of how these surfaces in 
configuration space look, in Fig. |T3. 

It is possible to define an ideal operator, Oideai, which will provide an ideal separatrix 
for distinguishing configurations which are in the domain of attraction of one or the other 
phase. Namely, we define 



C>ideai(config) 



traj . through config 



e(0^,(config(tmediurrr)) < (B)) . 



(A2) 



C^idcai of a configuration is the fraction of trajectories, starting at that configuration, which 
are in the metastable phase after time tmedium- For Hamiltonian dynamics, the measure on 



43 



the space of trajectories through a configuration is just the canonical measure of the tangent 
space; a point in the tangent space uniquely defines a trajectory. For Langevin dynamics, the 
measure for the space of trajectories through a point is given by the measure of realizations 
of the noise, with each trajectory corresponding to the noise realization which generates it. 
We do not use Cideai in our work because its measurement is impractical. 

The method presented in the main text for the determination of the nucleation rate, 
applied to an observable O, can be phrased as follows: first, we find the probability distri- 
bution as a function of the value of O; that is, we find the integral along each surface of 
constant O of the weight of the canonical ensemble. This gives 



p(i0-aHti<e/2) 

eP(0<ant) 



oc / d(surf. area) exp(-i//T) [rfO/rf(normal)]- , (A3) 

Jsurf 



the integral over the surface of the Boltzmann weight, divided by the surface normal deriva- 
tive of O. The surface normal derivative appears because, the larger the derivative is, the 
narrower is the region over which O differs by less than e/2 from Ocrit- 

Multiplying by {\AO/At\) precisely compensates for this normal derivative factor. This 
is because the mean value of \AO/At\, at any point in configuration space, is proportional 
to the gradient of O at that point. To see this, note first that only motion normal to 
the surface of constant O matters. Next note that, by our requirements on the dynamics, 
the rms. metric distance traveled in one direction, on averaging over possible momenta (or 
Langevin noise realizations), is independent of position in configuration space or direction 
considered. Therefore 



AC 



At 



^ /surf djsmi. area) exp{-H/T) 

Isuri ^isuri . area) exp(— iJ/T) [dO / d(jioTm.a\)]~ 



The product is proportional to the Boltzmann weighted area of the separatrix, and hence 
the fiux through the separatrix, as claimed in the body of the paper. 

The fiux through the separatrix will clearly differ for different choices oi O, as illustrated 



in Fig. [T^. It remains to show that d^, as defined in Eq. p.9| , precisely turns the correctly 
normalized fiux of trajectories through the separatrix into a count of trajectories which 
mediate nucleations. Fig. |l^ illustrates why this is true. First consider a trajectory which 
crosses the separatrix an even number of times and returns to the phase it came from, like 
Trajectory 1 in the figure. If the trajectory crosses n times, it gets counted n times in the 
sampling procedure, once at each crossing of the separatrix. Each count is with the same 
weight. For Hamiltonian dynamics this is because the evolution conserves energy (hence the 
Boltzmann factors are all the same) and phase space measure. For Langevin dynamics it 
is because the Langevin dynamics correctly generate the thermal ensemble. In either case, 
the number of times the trajectory gets sampled is in proportion to its contribution to the 
fiux, and each time it contributes zero to d; hence d correctly accounts for the number of 
nucleations (zero) the trajectory causes. 

Next consider a trajectory which does get from the metastable phase to the stable one. 
It is guaranteed to cross both separatrices an odd number of times, because each is to the 
right of B where the trajectory starts, and to the left of D where it ends. If the trajectory 
crosses a separatrix n times, it will appear in the average, used to determine d, n times, 
each with the same weight; and each appearance contributes 1/n to the determination of 
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^Trajectory 1 



"Trajectory 3 




Separatrix One 



"Trajectory 2 



Separatrix Two 

FIG. 16. cartoon of separatrices for two measurables, and some trajectories. Trajectory 1 
crosses separatrix 2 but not separatrix 1; but it (therefore) crosses an even number of times, and 
leads to a entry in determining dc)2- Similarly Trajectory 2 crosses separatrix 1 but not separatrix 
2, but also does not contribute to doi- Trajectory 3 does lead to a nucleation. It is sampled 3 times 
in computing do2, each time contributing (1/3); and sampled once in computing d^i, contributing 
1. Hence it gives the same contribution to the total rate computation for each measurable. 



d, so the number of nucleations is correctly counted as 1. Again, the Boltzmann weighted 
amount of flux the trajectory crossing represents, is the same at each of its crossings — 
of either separatrix. Hence if there is a larger total flux through, say, the separatrix for 
Oi, then this individual trajectory represents a smaller fraction of that flux, and gets an 
appropriately smaller weight in the sampling for determining doi. Its positive contribution 
to the value of doi is correspondingly smaller. Since the set of trajectories which mediate 
nucleations are the same, whichever separatrix we use to sample them, the smaller size of d 
exactly compensates for the larger flux through the separatrix. This is why the total rate 
computed is the same for each choice of O. 

While different observables give the same answer for the nucleation rate, they are not 
equally convenient as numerical tools. A good measurable must satisfy two conditions. First, 
it must be easy to measure it, so it can be used practically in reweighting configurations in 
the Monte Carlo (see subsection [IV B|) . Second, when determining d, statistics must accu- 
mulate with a reasonable sample of configurations, which means that a reasonable fraction 
of trajectories through the separatrix must actually lead to bubble nucleation. Otherwise it 
may take an exponentially large sample of trajectories to determine d with good statistical 
accuracy. Roughly, this will require that the separatrix of the observable is close to degener- 
ate with the separatrix of 0ideai- Note that, by construction, under Langevin dynamics half 
of trajectories through the Cideai separatrix lead to nucleations, and that this is the upper 
bound among all observables. 
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